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ABSTRACT 

We present a detailed dynamical analysis of the projected density and kinematical data available for 
the globular cluster lj Centauri. We solve the spherical anisotropic Jeans equation for a given density 
profile to predict the projected profiles of the RMS velocity a(R), in each of the three orthogonal 
coordinate directions (line of sight, proper motion radial, and proper motion tangential). The models 
allow for the presence of a central dark mass, such as a possible intermediate-mass black hole (IMBH). 
We fit the models to new HST star count and proper motion data near the cluster center presented 
in Paper I, combined with existing ground-based measurements at larger radii. The projected density 
profile is consistent with being flat near the center, with an upper limit 7 < 0.07 on the central 
logarithmic slope. The RMS proper motion profile is also consistent with being flat near the center. 
The velocity anisotropy profile, distance, and stellar mass-to-light ratio are all tightly constrained 
by the data, and found to be in good agreement with previous determinations by van de Ven et al. 
To fit the kinematics we consider anisotropic models with either a flat core (7 — 0) or a shallow 
cusp (7 = 0.05). Core models provide a good fit to the data with Mbh = 0; cusp models require 
a dark mass. If the dark mass in cusp models is an intermediate-mass black hole (IMBH), then 
Mbh = (8.7 ± 2.9) x 10 3 M Q ; if it is a dark cluster, then its extent must be < 0.16 pc. Isotropic 
models do not fit the observed proper motion anisotropy and yield spuriously high values for any 
central dark mass. These models do provide a good fit to the Gauss-Hermite moments of the observed 
proper motion distributions (/14 = —0.023 ± 0.004, he = 0.001 ± 0.004). There are no unusually 
fast-moving stars observed in the wings of the proper motion distribution, but we show that this does 
not strongly constrain the mass of any possible IMBH. The overall end-result of the modeling is an 
upper limit to the mass of any possible IMBH in uj Centauri: Mbh ^ 1-2 x 10 4 M© at ~ la confidence 
(or < 1.8 x 10 4 M at ~ 3c confidence). The la limit corresponds to M^ n /M tQt < 0.43%. We 
combine this with results for other clusters and discuss the implications for globular cluster IMBH 
demographics. Tighter limits will be needed to rule out or establish whether globular clusters follow 
the same black hole demographics correlations as galaxies. The arguments put forward by Noyola 
et al. to suspect an IMBH in lo Centauri are not confirmed by our study; the value of Mbh they 
suggested is firmly ruled out. 

Subject headings: globular clusters: individual (u> Centauri) — stars: kinematics. 



1. INTRODUCTION 

Intermediate-mass black holes (IMBHs) are of interest 
in a variety of astrophysical contexts (see, e.g., the re- 
views by van der Marel 2004; Miller & Colbert 2004). 
Especially the possible presence of IMBHs in the centers 
of globular clusters has intrigued astronomers for decades 
(e.g., Bahcall & Wolf 1976). This subject area under- 
went a considerable resurgence in recent years due to a 
combination of two advances. First, theoretical model- 
ing suggested plausible formation scenarios for IMBHs 
in the centers of globular clusters from realistic initial 
conditions (Portegies Zwart & McMillan 2002). Second, 
the quality of observational data sets, especially those 
obtained with the Hubble Space Telescope (HST), ad- 
vanced to the point where it is possible to measure ve- 
locity dispersions at a spatial resolution comparable to 
the sphere of influence size for plausible IMBH masses. 
This led to the identification of three globular clusters 
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with candidate IMBHs, namely M15, Gl, and u> Cen- 
tauri. 

Gerssen et al. (2002) used line-of-sight velocity data for 
individual stars in the core-collapsed cluster M15 (NGC 
7078) to argue for the presence of ~ 3000 M Q in dark ma- 
terial near the center, possibly in the form of an IMBH. 
Dull et al. (1997) and Baumgardt et al. (2003a) showed 
that this could be a plausible result of mass-segregation, 
and need not be an IMBH. Moreover, Baumgardt et 
al. (2005) argued that the observed steep brightness pro- 
file of M15, typical of the class of core-collapsed clusters, 
is in fact not consistent with expectation for a relaxed 
distribution of stars around an IMBH. Ho, Terashima, & 
Okajima (2003) and Bash et al. (2008) failed to detect 
X-ray or radio emission coincident with the cluster cen- 
ter. This is not necessarily inconsistent with an IMBH, 
since the accretion rate and/or the accretion efficiency 
are likely to be low. Nonetheless, M15 does not appear 
to be as strong an IMBH candidate as once thought. 

Gebhardt et al. (2002, 2005) used the velocity dis- 
persion measured from integrated light near the center 
of the M31 cluster Gl to argue for the presence of a 
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(1.8 ± 0.5) x 10 Mq dark mass at the cluster center. 
The relaxation time of Gl is too long to have produced 
such a concentration of dark mass due to mass segre- 
gation. Therefore, this result was taken as evidence of 
for the presence of an IMBH. Baumgardt et al. (2003b) 
disagreed with this interpretation, and argued that mod- 
els with no dark mass (IMBH or otherwise) provided an 
equally acceptable ht. The crux of this disagreement 
lies in the extremely small size of the sphere of influence 
(GMbh/c 2 ~ 0.03") of the putative IMBH, which is less 
than the spatial resolution of the observations. Phrased 
differently, the IMBH induces only a small signature in 
the data. As a result, the interpretation of the data is 
sensitive to the exact details of the data-model compari- 
son. The statistical significance of the IMBH detection is 
not high, especially if one takes into account the possibil- 
ity of "black hole bias" (discussed later in Section 16.3ft . 
On the other hand, the possible presence of an IMBH in 
Gl recently gained further credence with the detection 
of weak X-ray and radio emission from near the cluster 
center (Pooley & Rappaport 2006; Kong 2007; Ulves- 
tad, Greene & Ho 2007). The emission properties are 
consistent with accretion onto an IMBH, but other ex- 
planations cannot be ruled out. 

Noyola et al. (2008; hereafter NGB08) argued for the 
presence of an IMBH in the center of the cluster ui Cen- 
tauri (NGC 5139), which is well known for its many enig- 
matic properties (e.g., Meylan 2002). The dynamics of 
this cluster were studied previously by van de Ven et 
al. (2006; hereafter vdV06) using sophisticated axisym- 
metric models. However, there was insufficient kinemat- 
ical data at small radii to meaningfully constrain the 
possible presence of any dark mass. NGB08 performed 
integrated-light measurements of the line-of-sight veloc- 
ity dispersion for two 5 x 5" fields, one on the cluster 
center that they had determined, and one at R = 14" 
from that center. They obtained <7i os — 23. 0± 2.0kms _1 
and CTios = 18.6 ± 1.6kms _1 for these fields, respectively. 
The increase in cri os towards the cluster center was in- 
terpreted to imply the presence of an IMBH of mass 
M B h = 4.0t?;o M Q- The y also determined the sur- 
face brightness profile of lu Ccntauri from HST images, 
and inferred a shallow central cusp of logarithmic slope 
7 = 0.08 ± 0.03. This is consistent with theoretical pre- 
dictions for the cusp induced by an IMBH (Baumgardt 
et al. 2005). 

Both Gl and to Centauri are somewhat atypical globu- 
lar clusters. They are the most massive clusters of their 
parent galaxies, M31 and the Milky Way, respectively. 
Both clusters have been suggested to be the stripped nu- 
clei of larger galaxies (e.g., Freeman 1993; Meylan et 
al. 2001). Indeed, the nuclei of galaxies have very simi- 
lar structural properties to globular clusters (these nuclei 
are therefore also referred to as nuclear star clusters; e.g., 
Walcher et al. 2005). Accretion activity is known to exist 
in at least some of these nuclei (e.g., Seth et al. 2008). 
So if Gl and u> Centauri contain IMBHs, then these may 
simply be the low-mass tail of the distribution of super- 
massive black holes known to exist in galactic nuclei. As 
such, they would not necessarily be telling us anything 
fundamentally new about the formation and evolution 
of globular clusters in general. Nonetheless, Gl and u> 
Centauri are probably the best current candidates for 
globular clusters with IMBHs. Validating the existence 



of their suggested IMBHs is therefore of fundamental im- 
portance. 

The (initial) IMBH evidence for M15, Gl, and u Cen- 
tauri was all based on studies of line-of-sight velocities. 
Such studies are complicated because they require spec- 
troscopy in regions of high stellar density. Integrated- 
light measurements work well for a cluster as distant as 
Gl (Gebhardt et al. 2002, 2005), but the downside is 
that the sphere of influence in such clusters is poorly re- 
solved. In closer Milky Way clusters, the limited number 
of stars (especially brighter ones) introduces significant 
shot noise in integrated-light measurements. For uj Ccn- 
tauri, this required special treatment of the data to coun- 
teract this (NGB08). For M15, HST/STIS was used to 
obtain spectroscopy of individual stars (van der Marel et 
al. 2002). However, the signal-to-noise requirements for a 
velocity measurement limit the data set to the brightest 
stars. In turn, this sets the smallest radius inside which 
a velocity dispersion can be measured. 

A significant improvement in data quality is possi- 
ble with proper motion measurements. Such studies do 
not require much observing time and allow measurement 
even of the more numerous faint stars. Moreover, the 
ratio of the velocity dispersions in the tangential and 
radial proper motion directions directly constrains the 
velocity-dispersion anisotropy of the system. This is im- 
portant because of the well-known mass-anisotropy de- 
generacy for stellar systems (Binney & Mamon 1982). 
To determine accurate proper motions close to the cen- 
ter requires an observing platform with high spatial res- 
olution and long-term stability. The HST provides such 
a platform. McNamara, Harrison, & Anderson (2003) 
measured proper motions for Ml 5, which were modeled 
by van den Bosch et al. (2006) and Chakrabarty (2006). 
This confirmed the presence of dark mass near the cluster 
center, but did not shed new light on the issue of whether 
this is due to mass segregation or an IMBH. McLaughlin 
et al. (2006) obtained and modeled proper motions for 
47 Tuc, and used these to set an upper limit of 1000- 
1500 Mq on the mass of any IMBH in that cluster. 

In Anderson & van der Marel (2009; hereafter Paper 
I) we presented the results of a new observational HST 
study of lo Centauri. For the central few arcmin of the 
cluster we determined the projected number density dis- 
tribution of some 10 6 stars, as well as accurate proper 
motions of some 10 5 stars. The exact position of the 
cluster center was determined in three independent ways. 
First, we determined the center of the number density 
distribution using various methods, taking care to ac- 
curately account for incompleteness. Second, we deter- 
mined the kinematical center, being the symmetry point 
of the RMS proper motion velocity field. And third, we 
determined the center of unresolved light seen in 2MASS 
data. All determinations agree to within 1-2", yielding 
a final estimate with an uncertainty of ~ 1". 

The NGB08 center position, and hence the position 
of their central spectroscopic field, is at R = 12" from 
the center determined in Paper I (and coincidentally, the 
second off-center spectroscopic field they studied is also 
at R = 12"). The differences between the centers de- 
termined in these two studies is likely due to a combi- 
nation of subtle effects that biased the NGB08 analysis, 
as discussed in Paper I. As a result of this, the density 
cusp found by NGB08 was not measured around the ac- 
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tual cluster center. Also, we showed in Paper I that 
their use of unresolved light (as opposed to star counts) 
increased the impact of shot noise. In Paper I we de- 
termined the number density profile around the actual 
center of lo Centauri. It was found to be well matched 
by a standard King model, even though a shallow cusp 
may not be ruled out. 

In Paper I we also determined the RMS proper mo- 
tions of stars detected in the same 5" x 5" fields stud- 
ied spectroscopically by NGB08. For the canonical dis- 
tance D = 4.8 kpc (vdV06), the results translate to 
one-dimensional values tr pm = 18.9 ± 1.3 kms -1 and 
ff pm = 19.0 ± 1.6 km s -1 , for their "central" and "off- 
center" fields respectively. So we did not confirm the 
NGB08 finding that one of their fields has a higher dis- 
persion. The observation that there is no significant 
difference in proper motion dispersion between the two 
fields is consistent with our finding that they are both at 
the same radius R = 12" from the actual cluster center. 

Our results of Paper I did not confirm the arguments 
put forward by NGB08 to suspect the presence of an 
IMBH in lo Centauri. However, this does not mean that 
such an IMBH may not be present after all. The size and 
quality of our proper motion data set are superb. This 
provides the opportunity to study the central dynamics 
of lo Centauri at a level of detail that is unmatched by 
other clusters. We therefore present here a detailed study 
of the dynamics of lo Centauri, with the primary goal 
of exploiting the new data from Paper I to constrain 
the mass of any possible IMBH. We also include existing 
ground-based data in our data-model comparisons, to 
obtain the best possible results. 

We restrict our modeling efforts in the present paper 
to spherical models. This is a simplification, given that 
lo Centauri is elongated and rotating (in fact, unusually 
so for a globular cluster). However, both the elonga- 
tion (Geyer, Nelles, & Hopp 1983) and the rotation rate 
(vdV06) decrease towards the cluster center. Therefore, 
the assumption of sphericity is likely to be reasonable 
in the central region most relevant for constraining an 
IMBH. We find support for this throughout our paper by 
showing that spherical models yield an anisotropy pro- 
file, distance, and mass-to-light ratio that are in excellent 
agreement with the values obtained by vdV06. More so- 
phisticated models, such as those of the type described 
by vdV06, van den Bosch et al. (2006) and Chaname et 
al. (2008), have the potential to constrain other aspects 
of the structure of lo Centauri, such as its inclination, 
rotation properties, and the presence of any kinematic 
subcomponents. However, spherical models are sufficient 
to obtain robust constraints on the central mass distri- 
bution. 

The structure of this paper is as follows. Section [2] dis- 
cusses the general modeling approach. Section[3]presents 
the data on star count and surface brightness profiles 
used for the analysis. Section [5] discusses parameterized 
fits to these data. Section[5]presents the kinematical data 
used for the analysis. Section [6] presents the kinemati- 
cal data-model comparison, with a determination of the 
best-fitting model parameters. Section 17: discusses how 
the new IMBH results for lo Centauri fit in with our gen- 
eral understanding of IMBH demographics in globular 
clusters. Section [3] presents and discusses the conclu- 



sions. 

2. MODELING APPROACH 

2.1. Mass Density 

Imaging observations of a cluster yield an estimate of 
A t obs,x(-R), the radial profile of the surface brightness in 
a photometric band X, averaged over concentric circles 
of radius R around the cluster center. This quantity is 
related to the projected intensity / b s (i?) according to 

/zx[magarcsec~ 2 ] = -2.51og/ob s [ioPC~ 2 ]+21.572+M x, 

where Mq^x is the absolute magnitude of the sun in the 
given photometric band. Throughout this paper, units in 
which quantities are expressed are given in square brack- 
ets where relevant. In the presence of Ax magnitudes of 
foreground extinction, the extinction-corrected profile is 
given by 

I(R) = 10° AA *I ohs (R). (2) 

For a spherical cluster, I(R) is the line-of-sight integra- 
tion over the three-dimensional luminosity density j(r). 
This integration can be inverted through an Abel trans- 
form (e.g., Binney & Tremaine 1987) to yield 



i r°° 


dl' 


7T Jr 


dR 



(R) 



dR 



(3) 



To express this quantity in units of Lq pc -3 one needs 
to specify the cluster distance D. Sizes in arcsec and pc 
are related according to 



i?[pc] = ^[arcsec] 7rD[kpc]/648. 



(4) 



The mass density of stellar objects and their remnants 



is 



p(r) = [M/L}(r)j(r), 



(5) 



where \M/L](r) is the average mass-to-light ratio at each 
radius. Gradients in [M/L](r) can be thought of as aris- 
ing from two separate ingredients: (i) mass segregation 
among luminous stars; and (ii) mass segregation of the 
heavy compact stellar remnants (heavy white dwarfs, 
neutron stars, and stellar-mass black holes) with respect 
to the luminous stars. The first ingredient generally 
causes a shallow gradient in \M/L](r) over the scale of 
the entire cluster (e.g., Gill et al. 2008). This has little 
effect on the model predictions in the area surrounding 
the very center, where we wish to search for the sig- 
nature of an IMBH. By contrast, the second ingredient 
can in some clusters yield a strong increase in [M/L](r) 
near the very center (e.g., Dull et al. 1997; Baumgardt et 
al. 2003a) that mimics the presence of a dark sub-cluster. 

Our modeling approach allows for exploration of mod- 
els with arbitrary \M/L](r). However, we have found 
it useful and sufficient in the present context to restrict 
attention to a parameterized subset. Based on the pre- 
ceding considerations, we write the mass density as 



p(r) = p(r) + p da 



K r ) = Tj(r). 



(6) 



Here p(r) is the mass density obtained under the assump- 
tion of a constant mass-to- light ratio T. The excess pdark 
is treated mathematically as a separate component. For 
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Pdark we have found it convenient to use the parameter- 
ization 



The function f3(r) = 1 



J v 2 measures the anisotropy 



Pdark = (3M dark al rk /4^)/(r 2 + al rk ) 5/i , (7) 

corresponding to a Plummer density distribution. Here 
-ft^dark is the total excess mass in dark objects and Odark 
is the characteristic scale length of their spatial distribu- 
tion. The quantities T, Mdark and adark provide a con- 
venient way for exploring a range of plausible [M/L](r) 
functions using only three parameters. 

By assuming that T is independent of r, we neglect 
mass segregation among the luminous stars. For uj 
Centauri there are in fact several pieces of evidence 
that it has not yet undergone complete mass segrega- 
tion, consistent with its long half-mass relaxation time 
of ~ io9-96±o.03 yearg (McLaughlin & van der Marel 
2005). First, Merritt, Meylan & Mayor (1997) con- 
structed non-parametric dynamical models for uj Cen- 
tauri. They found that at the large radii sampled by 
their data, p(r) cx j(r) to within the error bars. Second, 
vdV06 constructed more general dynamical models for 
a wider range of data, and they too found that at all 
radii sampled by their study [M/L]0r) is consistent with 
being constant. Third, Anderson (2002) studied luminos- 
ity functions (LFs) at different radii in lu Centauri using 
data obtained with HST. He found that variations in the 
LFs are less than would be expected if mass segregation 
had proceeded to the point of energy cquipartition. And 
fourth, in Paper I we studied the proper motion disper- 
sion of stars along the main sequence. This showed that 
lower mass stars do move faster, but by a smaller amount 
than would be expected for energy equipartition. 

2.2. Gravitational Potential 

The total gravitational potential of the system is 

*tot = *(r) + $dark(r) - GMbhA, (8) 

where G is the gravitational constant and Mbh is the 
mass of any central IMBH. The gravitational potential 
corresponding to pOr) is 

p(r')r' 2 dr' + / p(r')r' dr' 



-4ttG 



(9) 



The Plummer potential corresponding to /?dark is 

$daxk(r) - -GA/ dark /(r 2 + a 2 ^ 2 . (10) 



The potential ^tot could in principle be expanded to 
also include an extended halo of particulate dark matter, 
but that is more relevant for galaxies. There is no evi- 
dence that star clusters are embedded in dark halos. For 
the specific case of ui Centauri, vdV06 showed that there 
is no need for dark matter out to at least two half-light 
radii. 

2.3. Intrinsic Velocity Moments 

The Jeans equation for hydrostatic equilibrium in a 
spherical stellar system is (e.g., Binney & Tremaine 1987) 



(11) 



dr r dr 

where vOr) is the number density of the stars under con 
sideration. Because of the spherical symmetry, v 2 . = in 



in the second velocity moments, where u 2 an = Vg = v^. 
Models with f3 = are isotropic, models with < 
are tangentially anisotropic and models with < f3 < 1 
are radially anisotropic. The models may in principle 
be rotating if v$ ^ 0. However, for the purposes of the 
present paper we do not need to address how the sec- 
ond azimuthal velocity moment = + v^ 2 splits into 
separate contributions from random motions and mean 
streaming. 

The measured kinematics in a cluster are generally rep- 
resentative for only a subset of the stars (e.g., stars in 
that magnitude range for which kinematics can be mea- 
sured). The number density of this subset is 



v(r)=j(r)f(r)/\(r), 



(12) 



where A(r) is the average luminosity of a star in the sub- 
set and fOr) is the fraction of the total light that the sub- 
set contributes at each radius. Compact stellar remnants 
play no role in this equation, since they contribute nei- 
ther to vOr) (their kinematics are not probed) nor to j(r) 
(they emit negligible amounts of light). However, fOr) 
and A(r) can vary slowly as a function of radius due to 
mass segregation among the luminous stars. The quan- 
tity A(r) may also be affected by observational selection 
effects that vary as a function of radius. Nonetheless, 
variations in A(r) will generally be smaller than those 
in f(r), as long as kinematics are measured for stars of 
similar brightness at each radius. Both fOr) and A(r) are 
only expected to have a shallow gradient over the scale 
of the entire cluster. These gradients are much smaller 
than those in jOr), which can vary by orders of magni- 
tude. Moreover, in the present paper we are mostly con- 
cerned with the small area surrounding the very center. 
We therefore assume that 

f(r)/X(r) = constant, (13) 

which implies that, as before, we neglect mass segrega- 
tion among the luminous stars. 

It follows from the preceding that one can substitute 
jOr) for v{r) in the Jeans equation (jlip . without actu- 
ally needing to know the value of the constant in equa- 
tion (fT3"|) . The Jeans equation is then a linear, first-order 
differential equation for jv% with solution 



3{r') 



dr 



„ / P(r") dr" , . 

Or') exp 2 / £ dr'. 

L J r r" i 

This is a special case of a more general solution for cer- 
tain classes of axisymmetric systems (Bacon, Simien & 
Monnet 1983). We consider anisotropy functions of the 
form 

P(r) = (P oo r 2 +p rl)/(r 2 +r 2 a ), (15) 

where r a > is a characteristic anisotropy radius, /3o < 1 
is the central value of /3(r), and Poo < 1 is the asymptotic 
value at large radii. If {3 — (3^ then f3(r) is constant. 
With this (30r), equation (fT4|) reduces to 



jv 2 (r) - 



j(r') 



dr 



(r') 



V 2 - 


i_ r 2" 
r ' a 


/3oc 


r ii 

rr \ 


j /r' 2 - 




r 2 H 


~ ^ a _ 




Vr^JI 


V r 2 n 


-r 2 J\ 



-,Po 



dr', (16) 
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which is always positive. 

The integral in equation IT6|) can be used with the 
equations given for j(r) and $tot(»") in Sections 12.11 
and 12.21 to allow numerical evaluation of v 2 (r) . This re- 
quires knowledge of the observed surface brightness pro- 
file Hobs,x(R) and the known extinction Ax- It also re- 
quires choices for the scalar model quantities D, T, Mbh, 
-Mdarkj Odark, Po, Poo and r a , which can all be chosen to 
optimize the fit to the available kinematical data. 

2.4. Projected Velocity Moments 

To calculate quantities that are projected along the 
line of sight, consider an arbitrary point P in the cluster. 
We adopt a Cartesian coordinate system (pmr, los, pmt) 
with three orthogonal axis and its origin at the cluster 
center C. The pmr- axis (proper motion radial) lies in the 
plane of the sky and points from the projected position 
of C to the projected position of P. The los-axis (line 
of sight) points from the observer to the point C. The 
pmt-axis (proper motion tangential) lies in the plane of 
the sky, and is orthogonal to the pmr- and los-axes in 
a right-handed sense. The cylindrical polar coordinate 
system associated with (pmr, los, pmt) is (R, 4>, z) and 
the spherical polar coordinate system is (r,8,<fi), with 
9 6 [— 7r/2, 7r/2], and 9 = n/2 corresponding to the pmt 
axis. The velocities in the (pmr, los, pmt) system satisfy 

Vpnu — v r cos <f> — Vct> sin <j>, 
v\ os = v r sin (f> + cos <fi, 

V pm t = V e . (17) 

Since vpu^ = in a spherical system, the second moments 
are 



"pmr {r, 4>) =v*(r)[l- p{r) sin 2 <p] , 
v L( r > <t>) =«£(»•) P- - P(r) c °s 2 4>], 



v 2 pmt (r,c/>)=v^r)[l-p(r)], 

where we have used that v 2 ^ = v 2 e = v 2 [l — P(r)]. 
The projected second moments satisfy 



(18) 



i(v 2 pmi )(R)- 



I(vL)(R) = 



Hv 2 pint )(R) = 



°° 2rjv*(r)[l - P(r) + P(r)(R/r) 2 )} dr 
i Vr 2 - R 2 ' 

00 2rf%(r)[l- P(r)(R/r) 2 } dr 



R _ Vr 2 - R 2 
2rjv*(r)[l - P(r)} dr 



(19) 



Vr 2 - R 2 

The derivations of these equations use the facts that 
cos 2 4> = (R/r) 2 and dlos = d\/r 2 - R 2 — rdrj\Jr 2 — R 2 . 
The factor of two arises from the identical contributions 
from points with los < and los > 0. The angle brackets 
around the squared velocities on the left-hand side indi- 
cate the combined effect of averaging over velocity space 
and density- weighted averaging along the line of sight. 
Formally speaking, the weighting is with number density 
if one is modeling discrete measurements, and with lu- 
minosity density if one is modeling integrated light mea- 
surements. However, in the present paper we are assum- 
ing that v[f) oc j(r), so the two cases are equivalent. The 
integrals in equation (|19[) can be evaluated numerically 
once the Jeans equation has been solved. 



2.5. Projected Velocity Distributions 

The normalized projected velocity distributions 
C(v, R) at projected radius R are generally different for 
the los, pmr, and pmt directions. The Jeans equation 
suffices to calculate the second velocity moments of these 
distributions. To calculate the distributions themselves 
it is necessary to know the full phase-space distribu- 
tion function. This is not generally straightforward (or 
unique) for anisotropic models. However, for an isotropic 
model the distribution function depends only on the en- 
ergy, and it can be uniquely calculated for given j(r) 
(being proportional to both i/(r) and p(r) under the 
present assumptions) and $ tot (r), with the help of Ed- 
dington's equation (Binney & Tremaine 1987). Because 
of the isotropy, the distribution C- lso (v, R) is independent 
of the choice of the velocity direction (los, pmr, or pmt) 
and satisfies, in analogy with equation (fl"9]) . 



I(R)C iso (v,R) 



«(«) 2rj(r)C ioc (v,r) dr 



Here r n 



R Vr 2 - R 2 

(v) is defined to satisfy V&(r max ) ~ ' 



(20) 
2 with 



\P = — <&tot' The quantity Ci oc (v,r) is the local velocity 
distribution (before projection along the line of sight) at 
radius r, which can be written with the help of Edding- 
ton's equation as (van der Marel 1994a) 



j(r)&ac(v,r) 



1 







/o 


Id®) 



*(r) 



\v 2 - 4" 
(21) 



The integrals in equations (j2"0")) and (|2"Tj) can be 
used with the equations given for j(r) and ^totM in 
Sections 12.11 and 12.21 to allow numerical evaluation of 
Ci so (v, R). The result is a symmetric function of v [i.e., 
£i so (v,R) = Aso( — v , R)]. This remains true if the dis- 
tribution function is given a non-zero odd part in the 
angular momentum L z around some axis (correspond- 
ing to a rotating model), provided that Ci so (v, R) is then 
interpreted as the distribution averaged over a circle of 
radius R. 

2.6. Observed quantities 

Equations (jT9|) and (|20| define velocity moments and 
velocity distributions at a fixed position on the projected 
plane of the sky. However, to determine these quantities 
observationally it is necessary to gather observations over 
a finite area. In integrated-light measurements this hap- 
pens naturally through the use of, e.g., apertures, fibers, 
or pixels along a long slit of finite width. In discrete mea- 
surements of individual stars this happens by binning the 
data during analysis. In either case, the corresponding 
model predictions are obtained as a density-weighted in- 
tegral over the corresponding area in the projected (x, y) 
plane of the sky. For example, the second line-of-sight 
velocity moment over an aperture is 

«>a P - f lK s )(R)dxdy/ f I{R)dxdy. (22) 



ap 



Similar equations apply to the other velocity moments 
and to the velocity distribution C\ so (v, R). In integrated 
light measurements one can also take the observational 
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point-spread-function (PSF) into account by convolving 
both I(x,y) and I(v 2 os )(x,y) with the PSF before evalu- 
ation of the integrals in equation (|22[) . 
In the following we write 

m = lv 2 W 2 a = (v 2 W 2 a <. = (v 2 W 2 
<J\os — V^los/ap ' °pmr — \" pmr /ap ) "pint — Wpmt/ap • 

(23) 

It is also useful to consider the quantity 

Spm = ((«pmr + «pmt)/2> a ^ = ^jlmr + *pmt)A (24) 

which corresponds to the one-dimensional RMS velocity 
averaged over all directions in the plane of the sky. Note 
that velocities in the plane of the sky are not directly 
accessible observationally. Instead of <7 pmr and <7 pmt we 
measure proper motions S pmr and £ pm t that satisfy 

^..[masyr- 1 ] = a... [ kms" 1 ] / (4.7404 D[ kpc]). (25) 

Transformation to velocities requires an assumption 
about the cluster distance D. 

The overlying bar in the preceding equations is used 
to differentiate the root-mean-square (RMS) velocity a 
from the velocity dispersion a. These quantities are equal 
only if the mean streaming is zero (which is always true 
in the pmr direction). 

2.7. Optimizing and Assessing the Fit to the Data 

In general, we will want to model line-of-sight data- 
points crios,i i Aai os ^ available for i = 1, . . . , N, as well 
as proper motion datapoints £ pm ...j ± A£ pm ... j avail- 
able for j — 1, ...,M (with each proper motion data 
point referring to either the pmr or the pmt direction). 
To calculate model predictions we use the known sur- 
face brightness profile fi bs,x(R) and the extinction Ax- 
The questions are then how to assess the fit to the data 
for given model parameters, and how to determine the 
parameters that provide adequate fits. 

2.7.1. Distance and Mass-to-Light Ratio 
Suppose that we have obtained model predictions <xi os ,i 
and Xpm...^ for the datapoints, using trial values D and 

T for the distance and mass-to-light ratio, respectively. 
The fit can then be improved by scaling of the model 
predictions, which corresponds to rescaling of the mass 
and distance scales of the model. Let Fi os and F pm be 
factors that scale the model predictions in the line-of- 
sight and proper motion directions, respectively. The 
best-fitting values of these scale factors and their error 
bars can be calculated using the statistics 

Xlo S = z2[(^o S ,i ~ F 1°sO'1o S , 1 )/Acti os , 1 ] 2 , 

i 

Xpm = /J[(^pm...,j - -FpmS pm ... ) j)/AS pm ...j] 2 . (26) 
3 

The best fit in each case is the value with the minimum 
X 2 (i.e., dx 2 /alF... — 0), and the 68.3% confidence region 
corresponds to A\ 2 = X 2 ~ Xmin — T The inferred scal- 
ings imply that the distance and mass-to-light ratio are 
optimized for D fit = DFi os /F pm and Tg t = TF los F pm . 
These fits apply at the given values of the other model 
parameters, some of which are also affected by the scal- 
ing. The model parameters Mbh and M da rk are pro- 
portional to F^ os /F pm , and the model parameter adark is 
proportional to Fi os /F pm . 



The error bars in Dfn and Tg t (at fixed values of the 
other parameters) follow from propagation of the uncer- 
tainties in .Fios and F pm . The errors in Fi os and F pm are 

uncorrelated, but the errors in Da t and T can be highly 
correlated or anti-correlated. We will not generally quote 
the correlation coefficient in this paper, but nonetheless, 
this should be kept in mind where necessary. 

The use of scalings as described here provides a means 
of reducing the parameter space of the models by two. 
Also, it directly yields the best-fit distance and mass- 
to-light ratio if the remaining parameters can be de- 
termined or eliminated independently. The parameters 
Mbh, -Mdark and adark can be determined as described in 
Section [2. 7. 21 They primarily affect the datapoints near 
the cluster center. So they can also be effectively elim- 
inated by restricting the definition of the x 2 quantities 
in equation (j26f to the datapoints at large radius. The 
anisotropy profile (3(r) can be determined as described 
in Section [2231 

2.7.2. Dark Mass 

Any dark component in our models can be either a 
cluster with parameters (Md ar k, adark) or an IMBH of 
mass -Mbh- The latter is mathematically (but not as- 
trophysically) equivalent to a cluster with parameters 
(Mdark = Mbh, adark = 0). It is therefore not possi- 
ble to discriminate between an IMBH and a dark cluster 
of sufficiently small size (within a resolution determined 
by the properties of the data) . Conversely, we only need 
to explore the (Mdark, adark) parameter space, and this 
will automatically tell us (at adark — 0) what the limits 
on an IMBH are. In the following we will use the term 
"dark mass" to refer generically to a dark component 
that could be either a dark cluster as in equation ([7]) or 
an IMBH. We have not explored models that have both 
a dark cluster and an IMBH, although that would not 
have been difficult. 

At given /3(r), we start by calculating initial guesses 

for D and T by fitting only the data at large radii. We 
then calculate for a range of trial values (Mdark, adark) 
the model predictions <7i s,i and £ pm ...,j for all the dat- 
apoints. For each pair of trial values we calculate Dfn 
and Tfit as in Section 12.7.11 by using the quantities x 2 os 
and Xpm> including observations at all distances from 
the cluster center. We then map the contours of x 2 = 
Xiob + Xpm in tnc two-dimensional (M dark , a da rk) param- 
eter space. The best fit (Mdark, fit, adark, fit) is obtained as 
the point of minimum x 2 ■ The uncertainties of the fit 
follow from standard A% 2 statistics, with the 68.3% con- 
fidence region corresponding to A% 2 = X 2 ~ Xmin < 2.3 
(Press et al. 1992). If there are combinations with 
Miark = within the confidence region, then a dark mass 
is not required to fit the data, for the assumed (3(r). If 
there are combinations with Mdark 7^ and adark = 
within the confidence region, then the presence of an 
IMBH is consistent with the data, for the assumed f3{r). 
If we assume that the dark mass is an IMBH, then its 
68.3% confidence region is given by A% 2 < 1 at adark = 0. 

The presence of an IMBH not only increases the RMS 
velocity towards the center, but it also changes the shape 
of the velocity distribution. The distribution measured 
through an aperture around the center will have power- 
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law wings that are significantly broader than those of a 
Gaussian distribution (van der Marel 1994b). This leads 
to the possibility of detecting very high velocity stars 
(Drukier & Bailyn 2003). It is therefore useful to deter- 
mine and analyze the observed velocity distribution of 
the stars around the center, to further assess the plau- 
sibility of the presence of an IMBH. The observed dis- 
tribution can be compared to the predicted distribution 
Aso,ap( w ) that can be calculated for an isotropic model. 
It is also useful to calculate the Gauss-Hermite moments 
(van der Marel & Franx 1993) of the observed velocity 
distribution as a function of distance from the center. A 
sufficiently massive IMBH causes an increase in hi to- 
wards the center (van der Marel 1994b). The absence of 
broad wings or a central increase in /14 can be used to 
obtain an upper limit to the mass of any IMBH. 

2.7.3. Amsotropy Profile 

Binney & Mamon (1982) showed that there are many 
different combinations of f3(r) and the enclosed mass 
M{< r) in a spherical system that can produce the 
same observed projected RMS line-of-sight velocity pro- 
file <j\ os {R). The key to constraining the mass distri- 
bution, and hence the possible presence of an IMBH, is 
therefore to constrain the anisotropy profile /3(r). This 
is not possible when only RMS line-of-sight velocities are 
available. However, (3(r) can be constrained when infor- 
mation is available on deviations of the line-of-sight ve- 
locity distributions (LOSVDs) from a Gaussian, as func- 
tion of projected position on the sky (van der Marel & 
Franx 1993; Gerhard 1993). This information is fairly in- 
direct, and to exploit it one must construct complicated 
numerical models that retrieve the full phase-space dis- 
tribution function of the stars (van der Marel et al. 1998; 
Cretton et al. 1999; Gebhardt et al. 2000; Valluri et 
al. 2004; Thomas et al. 2004; van den Bosch et al. 2008). 
This has been explored with considerable success in the 
context of studies of early-type galaxies (e.g., Cappellari 
et al. 2007). 

In the present context, high-quality proper motion 
data are available. It is then considerably more straight- 
forward to constrain f3(r). This was demonstrated by 
Leonard & Merritt (1989), who were the first to use 
proper motion data to examine in detail_the mass dis- 
tribution in a star cluster. The function [E pmt /E pmr ](i?) 
is a direct measure of anisotropy as seen in the plane of 
the sky. This ratio is independent of either the distance 
or the mass-to-light ratio scaling of the model. There is 
only a weak dependence on the presence of a dark mass 
(and this only affects the very central region). Therefore, 
the one-dimensional observed function [E pmt /E pmr ](i?) 
tightly constrains the one-dimensional intrinsic function 
/3(r). In fact, for a spherical system there is a unique re- 
lation between the two (Leonard & Merritt 1989). This 
is the primary reason why in the present context it is 
possible to obtain robust results based on something as 
simple as the spherical Jeans equation, while much more 
sophisticated techniques have become standard practice 
for galaxies. 

To build some intuitive understanding of the meaning 
of [E pmt /E pmr ](i?), consider the case in which /3(r) is in- 
dependent of radius and jv%(r) oc r~^. The relevant inte- 
grals in equation (TT9"]) can then be expressed analytically 
in terms of beta-functions. Their ratio is independent of 



radius R and simplifies to 

[E pmt /E pmr ] = £(1 - /?)/(£ - (3). (27) 

For small f3 (i.e., a system not too far from isotropy) we 
can use a first order Taylor approximation, which yields 

pm (R) = 0(1 - |), (28) 

where we have defined the function 

(3 pm (R) = 1 - {[E pmt /E pmr ](i?)} 2 , (29) 

in analogy with the definition of /3(r). This implies that 
the deviation of [E pmt /E pmr ] from unity is a factor (1— j) 

of the deviation of [^tan/^r] 1 ^ 2 from unity. (This uses 
the Taylor approximation (1 — /3] 1 / 2 ~ 1 — (/3/2), and 
similarly for /3 pm .) Hence, [E pmt /E pmr (i?)] is a "diluted" 
measure of the anisotropy that is present intrinsically, 
with dilution factor (1 — |). 

The approximations made above cannot be used in the 
core of a cluster, where £ = 0. The assumption that the 
system remains scale free all along the line of sight is then 
not accurate. However, in the outer parts of a cluster the 
approximation is a reasonably accurate guide, with the 
value of £ depending on the exact cluster structure. For 
example, an isothermal sphere has £ = 2 at large radii, 
while a Plummer model has £ = 6. The corresponding 
dilution factors (1 — |) are 0.50 and 0.83, respectively. 
Real clusters typically fall between these extremes. More 
extensive analytical calculations of /3 pm (i?) are presented 
in, e.g., Leonard & Merritt (1989) and Wybo & Dejonghe 
(1995). Either way, analytical calculations are mostly 
useful here for illustration. In practice we calculate the 
integrals in equation (TT9"]) numerically. 

The extent to which the velocity distribution of the 
stars (in either the los, pmr, or pmt direction) differs 
from a Gaussian provides secondary information on the 
amount of velocity anisotropy in the system. The model- 
ing technique used here is insufficient to exploit this infor- 
mation, since we are not calculating actual phase-space 
distribution functions for anisotropic models (although 
this would in principle be possible, e.g., Dejonghe 1987; 
Gerhard 1993; van der Marel et al. 2000). Nonetheless, 
we can observationally characterize the Gauss-Hermite 
moments as function of projected distance from the clus- 
ter center. Also, we can calculate the predicted velocity 
distributions Ci so (v,R) for an isotropic model, and from 
this the corresponding Gauss-Hermite moments. Com- 
parison of the observed and predicted quantities provides 
an independent assessment of possible deviations from 
isotropy. 

2.8. Numerical implementation 

To apply the formalism described above we augmented 
a software implementation used previously in, e.g., van 
der Marel (1994a). The code uses a logarithmically 
spaced grid, from very small to very large radii (so that a 
negligible fraction of the cluster mass resides outside the 
grid). Relevant quantities are tabulated on this grid. In- 
tegrals are evaluated through Romberg quadrature (e.g., 
Press et al. 1992) to high numerical precision. Tabulated 
quantities are interpolated using cubic splines for use in 
subsequent integrations. This avoids the need for cal- 
culation of higher-dimensional integrals. The derivatives 
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dl/dR and dj/<N> in equations ([3]) and (f2"Tj) . respectively, 
are calculated using finite differences. The derivative 
d&tot/dr in equation (|16[) can be manipulated analyti- 
cally so that numerical differentiation is not needed. The 
accuracy of the code was verified by reproducing analyt- 
ically known results for special potential-density pairs 
(e.g., Plummer models). 

3. PROJECTED DENSITY PROFILE DATA 

The primary input for our models is the projected ra- 
dial density profile. To characterize this profile over the 
full range of radii spanned by to Centauri, we have com- 
bined the ground-based measurements at intermediate 
and large radii collected in Trager, King & Djorgovski 
(1995), with the new HST measurements at small radii 
presented in Paper I. Much of the measurements come in 
the form of projected star count profiles ^ pro j(-R). Our 
equations of Section [2] are mostly expressed in terms 
of the observed (or implied) surface brightness profile 
Mobs(-R)- The profiles v pm j(R) and /i bs(-R) are uniquely 
related to each other (see eq. [3D] below), given our as- 
sumption to neglect mass segregation among the lumi- 
nous stars. 

3.1. Ground-based Data 

Trager et al. compiled literature data from various 
sources. We rejected from their compilation the data 
points to which they did not assign the maximum 
weight (unity in their notation). This leaves only the 
most reliable literature sources for u> Centauri, namely 
the integrated-light photo-electric measurements of Gas- 
coigne & Burr (1956) and Da Costa (1979), and the num- 
ber density measurements of King et al. (1968). All these 
authors presented averages for concentric circular annuli. 
Da Costa also presented separate measurements using a 
drift scan technique, which were subsequently calibrated 
to estimate the profile along concentric circular annuli. 
King et al. presented results for two separate plates from 
Boyden Observatory. We follow Trager et al. by referring 
to the different data sets with the abbreviations GBANN, 
DANN, DSCAN, ADH-1792, and ADH-1846; the latter 
two abbreviations refer to the IDs of the photographic 
plates analyzed by King et al. For our analysis we used 
the values provided by Trager et al., who calibrated the 
published measurements to a F-band surface brightness 
scale with a common zeropoint. We assigned to each data 
point an error bar of 0.142 mag, based on the procedure 
and analysis of McLaughlin & van der Marel (2005; see 
their Table 6). This corresponds to the scatter in the 
data with respect to a smooth curve. 

3.2. HST Data 

Paper I presented number density measurements for 
concentric circular annuli based on HST data. Consider- 
able effort was spent to accurately determine the cluster 
center, and to correct the results for the effects of incom- 
pleteness. Error bars were determined from the Poisson 
statistics of the counted stars. Paper I discussed number 
density profiles for various magnitude bins, and showed 
that there is no strong dependence on magnitude. For 
the analysis presented here we have used the profile for all 
stars in the HST images with i?-band instrumental mag- 
nitude < — 10 (this is defined in Paper I and corresponds 



to stars measured with signal-to-noise ratio S/N > 100). 
The choice of magnitude range was motivated by the de- 
sire to maximize the number of stars (i.e., to minimize 
the error bars on the density profile) while maintaining 
a similar mean magnitude (~ —12) as for the sample of 
stars for which high-quality proper motions are available 
from Paper I. 

To extend the profile from the Trager et al. compila- 
tion to the smaller radii accessible with the HST data 
it is necessary to calibrate the projected number density 
^proj(-R) t° a projected surface brightness /i bs,v(-R)- It 
follows from equations (J), ©> US)- and |T3]) that 

^obs y [mag arcsec -2 ] = — 2.51og(f pro j[arcsec -2 ]) + Z, 

(30) 

where the zeropoint Z is equal to 

Z = 21.572+M Q y + 2.51og{(//A[L ])(7rL>[kpc]/648) 2 }, 

(31) 

where M Q y = 4.83 (Binney & Merrifield 1998). Al- 
though (//A[Lq]) can in principle be estimated, we have 
found it more reliable to calibrate Z so that the sur- 
face brightness profiles from Trager et al. and Paper I 
agree in their region of overlap. We have done this by 
including Z as a free parameter in the fits described be- 
low, which yields Z — 18.21. For a canonical distance 
D = 4.8 kpc (vdV06) this implies that (//A[L Q ]) = 1.0. 
With this Z there is excellent agreement between the 
different datasets, as shown in Figure [T^,. 

NGB08 studied the surface brightness profile of unre- 
solved light in to Centauri using HST images. In prin- 
ciple, the advantage of using unresolved light over star 
counts is that one can trace the cumulative effect of all 
the unresolved faint stars, rather than counting only the 
less numerous resolved ones. However, we showed in Pa- 
per I that even at HST resolution, the unresolved light is 
dominated by scattered light from the PSF wings of only 
the brightest stars. So in reality, the inherent statistics 
of unresolved light measurements are poorer than those 
of star count measurements. Moreover, we showed in Pa- 
per I that the combined photometric and kinematic cen- 
ter of lu Centauri is actually 12" away from the position 
where NGB08 identified it to be. There is no simple way 
to transform their measurements over concentric annuli 
to a different center position. For these reasons, we have 
not included the NGB08 surface brightness datapoints in 
the dataset used for the dynamical modeling. 

4. PROJECTED DENSITY PROFILE MODELS 

It is important for the modeling to have a smooth rep- 
resentation of the projected intensity profile I{R), since 
the derivative dl/dR is needed in equation for the 
calculation of j(r). For the dynamical modeling we have 
used a parametric approach, although this can in princi- 
ple be done non-parametrically (e.g., Gebhardt & Fischer 
1995). We have fit the data using a smooth function of 
the form 

I ohs (R)=I b 2^y<* (R/b)^ 

[i + (i?/k)T ( "~ 7)/Q [i + (R/c) 5 y (e ~ v)/ (z2) 

This functional form has no particular physical signifi- 
cance, but it is useful because it can fit a wide range 
of observed profiles. The profile has a power-law cusp 
/ oc i? -7 at radii R <C b. It then breaks to a logarithmic 
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Fig. 1. — (a; top panel) V-band surface brightness profile of id Centauri in mag arcs ec - 2 versus radius in arcsec. Data points are from 
various sources as described in the text. The curve is our best model fit of equation (f >2t . which has a central cusp slope 7 = 0.05. (b; 
bottom panel) Residuals of the fit in the top panel. The small solid dots (black in the on-line version) are the HST star count data from 
Paper L The bigger open symbols are from the compilation of ground-based data in Trager et al. (1995), as described in the text: GBANN, 
hexagons (green); DANN, circles (red); DSCAN, squares (magenta); ADH-1792, triangles (cyan); ADH-1846, crosses (blue). 



slope 77 at R~ b, and further to a logarithmic slope e at 
R ~ c > b. The parameters a and S measure the sharp- 
nesses of the breaks. The intensity if,, corresponding to 
a magnitude fib through equation ([T]) , sets the overall in- 
tensity scale. The profile is similar to a so-called "nuker" 
profile (Lauer et al. 1995), but with two breaks instead 
of just one. We will refer to it as a generalized nuker 
profile. 

The best fit to the dataset presented in Section [3] 
was determined using a Levenberg-Marquardt iteration 
scheme (Press et al. 1992). It has parameters fib = 17.74, 
b = 183.1", c = 2158", 7 = 0.05, rf = 2.37, e = 10.00 (the 
maximum value allowed during the iteration), a = 2.29, 
and S — 1.70. The fit and its residuals are shown in Fig- 
ure Q] The quality of the fit is excellent. For the subset 
of the data from Paper I, the \ 2 = 130.5 for 125 dat- 
apoints with an RMS residual of 0.04 mag. The total 
Xmin = 169.3 for all 159 datapoints, corresponding to 
151 degrees of freedom. This is acceptable at the ~ la 
level for a x 2 distribution. The residuals appear gener- 
ally consistent with random scatter in line with the error 



bars, with little evidence for systematic trends. 

The best fit has a central cusp slope of 7 = 0.05. To 
determine the statistical error on this value we also per- 
formed a range of fits in which 7 was fixed a priori, while 
all other parameters were varied to optimize the fit. The 
central structure of our best fit for 7 = is shown in Fig- 
ure [21 for comparison to our overall best fit. The curve of 
Ax 2 = X 2 — Xmin as function of 7 was used to derive the 
error bar on 7, which yields 7 = 0.05 ± 0.02. 

When using a generalized nuker fit, a constant surface 
brightness core is inconsistent with the data at ~ 2.5<t 
confidence. However, this parameterized approach may 
not yield an unbiased estimate of the asymptotic slope 
7 at the smallest radii. The data at all radii are given 
equal weight in the x 2 01 the fit, but it is obviously only 
the data near the very center that contain most of the 
relevant information on 7. Also, the two breaks at radii 
b and c in the generalized nuker profile fit correspond to 
the breaks associated with the traditional core and tidal 
radii of u) Centauri. The parameterization does not have 
sufficient flexibility to reproduce any possible change in 
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Fig. 2. — V-band surface brightness profile of u> Centauri in mag arcsec" 2 versus radius in arcsec. The data and the symbols are the 
same as in Figure [T^,, but this figure focuses more on the central region. The solid curve (black in the on-line version) is the overall best 
model fit of equation (132 ft . which has a central cusp slope 7 = 0.05. The long-dashed curve (blue) is the best fit when the cusp slope is fixed 
a priori to 7 = 0. The dashed (cyan) and dotted (magenta) curves are the Wilson and King models, respectively, presented in McLaughlin 
& van der Marel (2005). These were tailored to fit the compilation of ground-based data in Trager et al. (1995). 



profile slope well within the core radius (r = 141_{ 3 
arcsec; McLaughlin & van der Marel 2005). 

An alternative and simpler approach to determine 7 
is to fit a straight line to the (log R, log J) data at 
R < #max- For i? max = 15", this yields 7 = 0.00 ± 0.07. 
The values of 7 inferred for other values of R max in the 
range 0-20" are in statistical agreement with this. This 
analysis indicates that the central data points are con- 
sistent with a flat core. 

An additional assessment of models with a flat core 
can be obtained from comparison of the new data with 
previously published models. McLaughlin & van der 
Marel (2005) presented detailed surface brightness fits 
for a large sample of galactic and extragalactic globu- 
lar clusters. They considered three types of models with 
decreasing amounts of diffuse outer structure, the power- 
law model, the Wilson (1975) model, and the King (1966) 
model. The power-law models are similar to the param- 
eterizations that we have used here, but with a more re- 
stricted subset of the parameters (namely, 7 = 0, a = 2, 
and c = 00 were kept fixed). The models from Wilson 
and King are both based on a parameterized isotropic 
phase-space distribution function. The Wilson prescrip- 
tion differs from the better known King models in that 
it uses a slightly different lowering of the exponential en- 
ergy distribution, so as to produce a model that has a 
more extended but still finite outer structure. McLaugh- 
lin & van der Marel discussed u> Centauri as a special 
example (their section 4.2.1) and found that it is well fit 
by a Wilson model, but less so by a King model. The 
central structure of both models is shown in Figure O 
The Wilson model is very similar to our best fit model 
with fixed 7 = 0. The King model is also similar at radii 
R < 100", but it under-predicts the brightness at inter- 
mediate radii (before improving again at larger radii, not 
shown in Figure [2]). 

The Wilson and King models shown in Figure [2] are 
the ones that were previously fit by McLaughlin & van 
der Marel to the Trager et al. data. At small radii these 



models represent an extrapolation of the ground-based 
data under the assumption of a constant surface bright- 
ness core. Although the models were not fit to the HST 
data, they do provide a perfect fit to the two innermost 
HST star count data points. By contrast, our own best 
generalized nuker model fit with 7 = 0.05 over-predicts 
both of the innermost datapoints by about la. The fact 
that 7 = 0.05 is preferred in our fits over 7 = is there- 
fore because it fits better the data at radii R » 10", and 
not because it provides the best fit closest to the center. 

NGB08 inferred from their HST study of integrated 
light in lu Centauri that it has a central cusp slope 7 = 
0.08 ±0.03. At first glance, this result appears consistent 
with that from our generalized nuker fit. However, since 
they did not measure the brightness profile around the 
same cluster center, the agreement between our inferred 
cusp slopes is likely coincidental. 

In summary, we conclude that the central cusp slope 
of co Centauri is 7 < 0.07 at ~ la confidence. That is, 
we interpret the 7 = 0.05 ± 0.02 from our generalized 
nuker fit as an upper limit. We do this because it over- 
predicts the central two data points, and because other 
statistics do not show evidence for a non-zero slope. In 
the dynamical analysis that follows we specifically con- 
sider two generalized nuker fits that span the range of 
possibilities: models with a shallow cusp (7 = 0.05) and 
models with a constant density core (7 = 0). We will 
refer to these as cusp models and core models, respec- 
tively. In all calculations we have transformed the ob- 
served intensity I bs,v(R) to an extinction-corrected in- 
tensity iy(-R) using equation ([2]) with Ay = 0.372 for lu 
Centauri (McLaughlin & van der Marel 2005). 

5. KINEMATICAL DATA 

To constrain the parameters of the dynamical models 
we need to use both line-of-sight and proper motion data 
at a range of radii. For this we have combined ground- 
based line-of-sight and proper motion data at interme- 
diate and large radii with the new HST proper motion 
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data at small radii presented in Paper I. 

We construct spherical models to interpret the data. 
These models yield only a single radial profile for each 
projected quantity, with no azimuthal variations from the 
observed major axis of uj Centauri to its observed minor 
axis. To get accurate results for the model parameters, 
it is important that these models be used to reproduce 
the average radial variations seen in the data. This is 
achieved by using for each projected quantity that is fit 
the average of the measurements along concentric annuli. 
This is indeed what we extracted for the projected den- 
sity in Section [U and we therefore do the same here for 
the projected velocity distributions. When averages are 
taken over annuli, all the odd moments of the velocity 
distribution are zero. Of primary interest here are the 
second moments, which yield the RMS velocity cti os and 
the RMS proper motions S pm r and S pm t- 

In reality, variations do exist between observed quan- 
tities along the major and minor axis. vdV06 discussed 
this in detail for the existing ground-based data. In Ap- 
pendix B we address this for the new HST proper motion 
data. In particular, the RMS proper motions are not the 
same along the major and minor axes, and the RMS 
proper motions in the major and minor axis directions 
on the sky are not the same. Also, the mean proper mo- 
tion in the transverse (pmt) direction is not zero (vdV06; 
Section f5. 2. 31 Appendix C). By modeling these issues it 
is possible to constrain, e.g., the cluster inclination and 
its rotation rate. However, it is not essential to model 
these issues to achieve a good understanding of the cen- 
tral mass distribution, which is the main goal here. 

5.1. Ground-based Data 
5.1.1. Line-of-Sight Kinematics 

vdV06 compiled line-of-sight velocity data of individ- 
ual stars in lo Centauri. The original data were obtained 
by Suntzeff & Kraft (1996), Mayor et al. (1997), Rei- 
jns et al. (2005), and Gebhardt et al. (in preparation). 
vdV06 performed various operations and cuts to homog- 
enize and correct the data as necessary for dynamical 
modeling. This included: improved removal of cluster 
non-members; removal of stars with low enough bright- 
ness or large enough velocity error bars to make their 
reliability suspect; application of additive velocity shifts 
to correct for offsets in velocity zero-points between data 
sets; multiplicative rescaling of the velocity error bars 
where necessary, as dictated by actual measurements of 
scatter between repeat measurements; subtraction of the 
systemic line-of-light velocity; and subtraction of the ap- 
parent solid-body rotation introduced into the velocity 
field through the combination of systemic transverse mo- 
tion and large angular extent ("perspective rotation"). 
This produced a final homogenized sample of line-of- 
sight velocities for 2163 stars with uncertainties below 
2.0 kms" 1 . 

Glenn van de Ven kindly made the homogenized sam- 
ple from vdV06 available to us for use in our data-model 
comparisons. To determine the profile of a\ as (R) we 
binned the stars in radius and then used the statistical 
methodology described in Appendix A. Radii were calcu- 
lated with respect to the center determined in Paper I, 
which differs by 15.5" from the center used by vdV06 
(which was determined by van Leeuwen et al. 2000). We 
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chose a central bin of 20" radius (with 25 stars) for max- 
imum resolution near the center, and then chose the re- 
maining bins to contain 200 stars each. 

Sollima et al. (2009) obtained line-of-sight velocities of 
318 cluster members at large radii, which they merged 
with observations closer to the center by Pancino et 
al. (2007). This yielded a homogeneous sample of 946 
cluster members, from which they derived the profile 
oios(-R)- Seitzer (1983) presented the CTi os (i?) profile de- 
rived from 118 clusters members with known velocities 
at a range of radii. We included the published profiles 
from both studies in our data-model comparisons. For 
each study we averaged the first two bins in the pro- 
file together, since the first bin had a large error bar and 
did not provide spatial resolution that competes with the 
vdV06 compilation. 

Scarpa et al. (2003) obtained line-of-sight velocities of 
75 cluster members at large radii. From this they de- 
termined the velocity dispersion cri os in four radial bins. 
The stars were limited to the West side of the cluster, 
so the results do not represent averages over a complete 
annulus on the sky. Nonetheless, we included these obser- 
vations in our data- model comparisons. We subtracted 
in quadrature from the published uios values the average 
measurement error of 1 kms -1 . To obtain a\ os we added 
in quadrature V 2 /2, where V is the rotation velocity am- 
plitude quoted in Scarpa et al. (2003); the factor 1/2 is 
the average of a squared sinusoidal variation over a full 
range of angles. 

The profiles obtained from the Sollima et al. (2009), 
Seitzer (1983), and Scarpa et al. (2003) observations do 
not pertain to the same center as that inferred in Paper I. 
However, the smallest radial bin that we use from these 
studies is the 4 arcmin diameter combined central bin of 
Sollima et al. A center offset of 15.5" (as applicable to 
the Sollima et al. and Scarpa et al. results, who used the 
van Leeuwen et al. 2000 center) should not significantly 
affect statistics calculated over bins that extend to such 
large distances from the center. 

At small radii we used as additional data points the 
integrated- light measurements from NGB08. They in- 
ferred velocity dispersions over square fields of 5" x 5". 
Because of the significant size of these fields, we did not 
include the impact of PSF convolution in our model pre- 
dictions for these data. We placed the measurements at 
the distance from the cluster center inferred in Paper I, 
which is 12" for both fields. We cannot use the data to 
calculate a true average over an annulus. However, we 
do have two measurements at the same distance from the 
center at a 60° angle. Including both measurements in 
our model fits is in low-order approximation similar to 
fitting the average over an annulus. 

NGB08 did not measure the mean velocity over their 
fields as compared to the systemic velocity. Hence, their 
measurements are dispersions cri os rather than RMS ve- 
locities trios- However, it is known from previous observa- 
tions (e.g., vdV06) that V\ os /cr\ oa in u Centauri decreases 
towards the center. This ratio is small enough at the 
position of the NGB08 fields that <7i os <7i os to high 
accuracy (see Section 15.2.31 for a quantitative estimate). 

Figure [3^ shows the collected <7i os data from all the 
sources, as function of projected distance R from the 
cluster center. The different data sets match well in their 
regions of overlap. 
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Fig. 3. — Observed kinematical quantities for ui Centauri as function of projected radius R from the cluster center. Data points are 
based on analysis of various sources described in the text: circles (magenta in the on-line version) are based on vdV06, which is itself a 
compilation of line-of-sight data from Suntzeff & Kraft (1996), Mayor et al. (1997), Reijns et al. (2005), and Gebhardt et al. (in prep.), and 
proper motion data from van Leeuwen et al. (2000); triangles (cyan) are based on NGB08; stars (blue) are based on Sollima et al. (2009); 
squares (green) are based on Seitzer (1983); crosses (red) are based on Scarpa et al. (2003); and small dots (black) are based on the 
HST data from Paper I. Some error bars are of similar size as the plot symbols. Data from all ground-based sources except NGB 08 wer e 
multiplied by a factor 1.023 to correct to same the characteristic main-sequence stellar mass as for the HST sample (see Section 15.1.31 . 
(a; top panel) RMS line-of-sight velocity <ti os in km/s (data indicated with open symbols), (b; middle panel) RMS proper motion S pm in 
mas/yr, averaged over all directions in the plane of the sky, as defined by equation 1241 (data indicated with closed symbols), (c; bottom 
panel) Ratio S pm t/S pmr of the RMS motions in the tangential and radial proper motion directions (data indicated with closed symbols). 
The RMS quantities in each dataset are averages over adjacent circular annuli. The extent of the annulus for data point i can be assessed 
by taking half the distance between adjacent datapoints, ARj PS (iJi+i — Ri-\)/2. 



5.1.2. Proper Motion Kinematics 

Ground-based proper motion data are available from 
van Leeuwen et al. (2000). vdV06 performed various op- 
erations and cuts to optimize and correct the data as 
necessary for dynamical modeling. This included: im- 
proved removal of cluster non-members; removal of stars 
with low enough brightness, large enough velocity er- 
ror bars, or neighbors that are nearby enough to make 
their reliability suspect; and subtraction of the apparent 
solid-body rotation introduced by perspective rotation, 
and other spurious solid body rotation components. This 
produced a final sample of proper motions for 2295 stars 
with uncertainties below 0.2masyr _1 (i.e., 4.6 km s -1 at 
the canonical distance D — 4.8 kpc). 

Glenn van de Ven kindly made the homogenized sam- 



ple from vdV06 available to us for use in our data-model 
comparisons. Proper motions were given along the major 
and minor axis directions. For each star we transformed 
this into proper motions in the radial (pmr) and tangen- 
tial^ (pmt) directions. We then determined the profiles 
of E pmr (i?) and S pm t(i?) by binning the stars in radius 
and using the statistical methodology described in Ap- 
pendix A. Radii and proper motion coordinate transfor- 
mations were calculated with respect to the center deter- 
mined in Paper I. We chose a central bin of 90" radius 
(with 30 stars) for maximum resolution near the center, 
and then chose the remaining bins to contain 200 stars 
each. We followed vdV06 by excluding the 65 outermost 
stars with radii R ~ 20-30 arcmin. 
For fitting models to the data we always use the quan- 
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tities E pmr and E pmt . However, for visualizing the results 
we have found it more convenient to use the quantities 
E pm and E pm t/E pmr . The RMS proper motion E pm aver- 
aged over all directions in the plane of the sky can be cal- 
culated directly through its definition in equations (124| 
and (|25p . In practice, we found it more convenient to 
calculate it by applying the statistical methodology de- 
scribed in Appendix A on the combined array of the 
pmr and pmt data points. For visualization of the ra- 
tio E pmt /E pmr we adopted a radial binning scheme with 
400 stars in each bin, to obtain smaller and more useful 
error bars. 

Figure [3b shows the E pm data as function of pro- 
jected distance f?_from the cluster center, and Figure [3t 
shows the E pmt /E pmr data. The velocity distribution is 
close to isotropic near the center, but with some radial 
anisotropy. There is a significant increase in tangential 
anisotropy at a few core radii. The latter is mostly due to 
the contribution of rotation to the second moment E pmt . 
The anisotropy in the proper motion dispersions is much 
less (King & Anderson 2002). 

5.1.3. Dependence of Kinematics on Stellar Mass 

Luminous stars in a cluster undergo two distinct effects 
as a result of two-body relaxation: mass segregation and 
energy equipartition. These effects are intimately con- 
nected through basic dynamical theory, but from an ob- 
servational perspective it is useful to think of them as 
separate effects. Mass segregation causes stars of differ- 
ent masses to have slightly different spatial distributions. 
In Section [5] we assumed that this could be neglected 
when parameterizing the run of mass-to-light ratio with 
radius, and when connecting star count data at small 
radii to surface brightness data at large radii. This as- 
sumption was motivated by the absence of high-quality 
data to constrain the amount of mass segregation over 
the full scale of the cluster, and by the desire to mini- 
mize the number of a priori model assumptions. Energy 
equipartition causes stars of low mass to move faster than 
those of high mass. We measured in Paper I that low- 
mass stars in to Centauri do indeed move faster than 
high-mass stars, albeit by less than would be predicted 
for complete energy equipartition (v oc ?7i~ - 5 ). Since 
we have a direct measurement of this, there is no reason 
to ignore this effect in the modeling. We therefore cor- 
rect for it when comparing the observed kinematics from 
different data sets. 

The available discrete ground-based line-of-sight and 
proper motion data all pertain to relatively bright stars. 
The apparent magnitudes of these stars place them on 
the giant or sub-giant branch. Since stellar evolution 
proceeds relatively rapidly after the main-sequence turn- 
off, these stars all have essentially the same mass as the 
main-sequence turn-off. By contrast, the stars in our new 
HST proper motion sample from Paper I (discussed in 
detail in Section [5. 21 below) extend from about the main- 
sequence turn-off to ~ 3 magnitudes below that. From 
the analysis of Paper I it follows that the RMS proper 
motion of stars at the average magnitude of the HST 
sample is ~ 2.3% higher than that at the main-sequence 
turn-off. 

For a fair comparison to the new HST proper motion 
data, we took the kinematics from the discrete velocity 
studies in Sections 15.1.11 and 15.1.21 and multiplied them 
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by a factor 1.023. This implies that we use the observed 
(sub-)giant star kinematics to estimate the kinematics 
of lower-mass stars in the same general population. We 
assume implicitly that the variation of RMS kinematics 
with mass is the same everywhere in the cluster (although 
it was only measured near the cluster center) , and that it 
applies also in the line-of-sight direction (although it was 
only measured in the proper motion direction). In real- 
ity, one would expect the relation between mass and ve- 
locity dispersion to change with radius, since relaxation 
proceeds faster near the cluster center than further out. 
This was neglected in the present context by applying 
the same small correction throughout the cluster. This 
was motivated in part by the fact that our study focuses 
mostly on the central part of the cluster anyway, but also 
by the fact that we didn't actually detect a difference in 
equipartition in Paper I between the different (central 
and major axis) fields for which we obtained proper mo- 
tions. We did not apply the multiplicative correction to 
the NGB08 data. Their measurements are based on in- 
tegrated light, excluding the brightest regions in their 
fields. It is therefore not obvious that their measure- 
ments would be dominated by (sub-)giant stars. Either 
way, a 2.3% correction would be much less than the error 
bars on their datapoints. 

In summary, we are using our Jeans models to fit the 
kinematics of stars with masses typical of our HST proper 
motion sample. We extend our actual HST kinemat- 
ics to larger radii using ground-based kinematics that 
have been corrected to the same characteristic mass. We 
model the cluster density profile using HST star counts 
for a sample that is also centered on the same charac- 
teristic mass (see Section l3~2"T) . We make the simplifying 
assumption that there is no mass segregation among the 
luminous stars only so that we can extend the HST den- 
sity profile to larger radii, where it can be tied to ground- 
based star-count data and integrated photometry that 
are based largely on giant stars (see Section 13. ip . Since 
this approximation does not affect the very central region 
of the cluster, it should have little effect on our results for 
the central mass distribution. In fact, the multiplicative 
correction of 1.023 discussed above also has little effect 
on this. We verified explicitly that none of the main con- 
clusions of our paper change at a level that exceeds the 
formal uncertainties when the multiplicative correction 
is omitted. 

5.2. HST Proper Motion Data 
5.2.1. Sample selection 

In Paper I we derived proper motions for stars ob- 
served in two fields observed by HST, a "central field" 
on the cluster center, and a "major axis field" posi- 
tioned adjacent to it in non-overlapping fashion roughly 
along the cluster major axis. Here we use the "high- 
quality" sample presented in Paper I, which contains the 
non-saturated stars that are isolated enough and bright 
enough (£?-band instrumental magnitude < —11) for a 
particularly reliable proper motion measurement. This 
sample has 53382 stars in the central field and 19593 stars 
in the major axis field. In this sample, 95.1% of the stars 
have proper motion errors below 0.2 masyr™ 1 . Our cut 
in proper motion accuracy is therefore not very differ- 
ent from that applied by vdV06. However, the median 
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proper motion error per coordinate in the final sample 
of vdV06 is ~ 0.14 masyr -1 , whereas it is a factor two 
smaller, 0.07 masyr" 1 , for our sample. We note that our 
use of only the high-quality sample is rather conservative. 
Even the remainder of the sample of Paper I is quite ac- 
curate, and could well have been used to constrain the 
cluster dynamics. 

The central HST field covers radii R < 147.4", while 
the major axis field covers radii 99.7" < R < 347.1". 
The central field has complete position angle coverage 
out to R = 71.7" (the radius of the largest enclosed cir- 
cle). For 71.7" < R < 147.4" the coverage over the range 
of position angles becomes increasingly sparse. The ma- 
jor axis field is restricted at all radii to position angles 
within 48.5° of the major axis. For the analysis of the 
present paper it is important to have coverage of all posi- 
tion angles at a given radius, so that average kinematical 
quantities over circular annuli can be calculated. For this 
reason we have used only the 25194 stars at R < 71.7" 
to constrain the dynamical models. 

In our analysis we ignore the 47781 stars with high- 
quality HST proper motions at radii 71.7" < R < 347.1". 
We do discuss relevant quantities derived from these data 
in Appendix C. More sophisticated axisymmetric model- 
ing techniques (such as those in vdV06; van den Bosch et 
al. 2006; and Chaname, Kleyna, & van der Marel 2008) 
that might be applied in the future probably would not 
want to exclude these data in their analysis. By contrast, 
to include these data in our spherical models would re- 
quire estimates of the ratios 1Z of the average of either 
Spmr or S pm t over a restricted range of position angles, 
to its average over a circle. We have experimented with 
such estimates, but concluded that the approximations 
that need to be made introduce sufficient uncertainties 
that these data do not really help to constrain our mod- 
els. Conversely, we have found no evidence that inclusion 
of the HST data at 71.7" < R < 347.1", combined with 
reasonable estimates for the ratios 1Z, would alter in any 
way the conclusions draw below (see Appendix C). 

To reject cluster non-members we created scatter plots 
of total two-dimensional proper motion \p\ versus radius 
R, where p is the proper motion vector of an individual 
star. Such a plot reveals a sparse population of field 
stars at large values of \p\, in addition to the numerous 
cluster stars at low \p\. We rejected those stars from the 
sample that reside at R > 10" and for which \p\ > 2.8 — 
(r/500") masyr" 1 . At the radii for which we have HST 
proper motions, this equation is a good approximation 
to the escape velocity curve for a spherical model with a 
canonical distance and mass-to-light ratio (D = 4.8 kpc, 
Ty = 2.5; vdV06). With this criterion, only 27 non- 
member stars were identified at R < 71.7", the closest 
of which resides at R = 21.4". Since the presence of an 
IMBH may result in the presence of rapidly moving stars 
at small radii, we did not reject any stars at R < 10". 
The wings of the observed velocity distribution at small 
radii are discussed in detail in Section [6. 71 b elow . 

Our final HST proper motion sample for the dynamical 
modeling consists of 25167 stars. This is a factor ~ 11 
larger than the number of stars with proper motions in 
the vdV06 sample. Most of their stars reside at larger 
radii than those sampled here. At the radii R < 71.7" of 
our sample, vdV06 had only some two dozen stars with 
proper motions, and ~ 270 stars with line-of-sight veloc- 



ities. The HST data therefore provide a major advance 
for constraining dynamical models, especially as it per- 
tains to the gravitational influence of a possible IMBH 
at small radii. 

5.2.2. RMS Proper Motions 

For each star in the HST sample we decomposed the 
observed proper motion vector p into components in the 
pmr and pmt directions, respectively. We then binned 
the stars in radius and derived the profiles £ pmr (-R), 
^pmt(-R)) and E pm (i?). We did this in similar fashion 
as for the ground-based proper motion data discussed in 
Section l5.1.2[ using the statistical methodology described 
in Appendix A. 

Figure [5b shows the £ pm data as function of projected 
distance R from the cluster center. The results match 
well with those from vdV06 in the region where the_data 
sets overlap. For the individual data points of S pmt , 
£ pmr and E pm we used binning in annuli that are 1" 
wide. For the central aperture we used a circle of ra- 
dius 1.55", which gives an average radius R = 1.1" for 
the datapoints enclosed by it. This central aperture is 
larger than the ~ 1" uncertainty in the position of the 
cluster center derived in Paper I. This uncertainty should 
therefore not significantly affect the inferred kinematics. 
The adopted binning gives acceptably small error bars in 
all apertures, without losing too much spatial resolution 
near the center. 

One immediate result from the inferred S pm (R) profile 
is that there is no strong increase towards the center, con- 
sistent with arguments presented in Paper I. A straight 
line fit to the (R, S pm ) data for R < i? m ax, with i? max = 
15", yields a slope of -0.06 ± 0.08 km s" 1 arcsec" 1 . The 
slopes inferred for other values of i? m ax in the range 0- 
20" are statistically consistent with this. This analysis 
indicates that the central data points are consistent with 
a flat profile of RMS velocity as function of radius. 

The error bars on the observed RMS velocities can be 
reduced by increasing the amount of binning. For the 
ratio of S pmt /E pmr shown in Figure [3J: we used binning 
in annuli of 8" wide, so that even small deviations from 
isotropy can be measured. As for E pm , the results match 
well with those from vdV06 in the region where the data 
sets overlap. The average over all HST data points is 
S pmt /S pmr = 0.983±0.006. Hence, there is a very small, 
but significant, amount of radial anisotropy in the central 
arcminute of lo Centauri. 

5.2.3. Proper Motion Rotation 

Our data calibration of Paper I (see discussion in Sec- 
tion 3.6.4 of that paper) removed by necessity any possi- 
ble solid-body rotation component from the data. This 
is because we did not have access to stationary back- 
ground sources in the HST fields, as would be required 
to measure absolute motions. This calibration has the 
advantage that it automatically removes any perspective 
rotation present in the actual proper motions (however, 
this is negligible for the small radii R < 71.7" sampled 
by our data anyway; < 0.0035 masyr" 1 = 0.08 kms -1 
at the canonical D = 4.8 kpc). But of course, it has the 
disadvantage that it also removes a solid-body rotation 
fit to the actual rotation field of the cluster. This is a 
fundamental limitation of the data. However, we discuss 
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below that this has only negligible importance for our 
models. 

In Paper I, we also used a "local correction" that re- 
moved all other mean systematic motions (as opposed to 
random motions) from the proper motion data. However, 
that was by choice, and not by necessity. The motiva- 
tion for this is to calibrate out small time-variations in 
the higher-order geometric distortion terms. This im- 
proves the measurement of the random motions in the 
cluster. However, one would of course not want this to 
remove an important signal that is present in the data, 
and in particular, any possible differential rotation in the 
plane of the sky (i.e., the part of the rotation field that 
remains after subtraction of a solid-body fit). We tested 
for such differential rotation in Paper I, and found that 
it is negligible in the central HST field. Proper motion 
rotation curves (the mean tangential motion as function 
of radius) derived from data with and without the local 
correction agree to within 1 kms" 1 at all radii. There- 
fore, solid-body rotation is the only unrepresented rota- 
tion component that could affect the modeling. 

In principle, the loss of information on solid-body ro- 
tation can be partly recovered as in vdV06. They use 
the fact that for an equilibrium axisymmetric system the 
equation 

Vi os = 4.7404 D[ kpc] tani V y (33) 

must hold everywhere on the projected plane of the sky. 
Here i is the inclination, y the projected minor axis di- 
rection, and the symbol V is used to distinguish a mean 
proper motion from a mean velocity V. However, to use 
this equation with accuracy for our HST proper motions 
would require many more stars with line-of-sight veloci- 
ties than are actually available at R < 71.7". 

The pmr direction is perpendicular to the direction of 
absolute rotation, and our estimates of £ pmr are there- 
fore not affected by our inability to measure absolute 
rotation. Also, in an equilibrium system there cannot 
be mean streaming in the pmr direction, when aver- 
aged along a circle. However, our calculated values of 
Sp m t = (V pmt + £p mt ) 1/2 wil1 be systematically low by 
a multiplicative factor g < 1, because we are unable to 
include any actual solid-body contribution from V pmt . 
The proper motion data presented by vdV06 did include 
measurements of absolute rotations in the plane of the 
sky. Figure [4] shows for their data the ratios V pm t/£ pm t 
and g = E pmt /£ pm t, where the quantities in the nomi- 
nator and denominator are all averages over an annulus. 
These measurements allow us to estimate the value of 
g for the HST data. The central data point shows ex- 
plicitly that there is very little rotation in w Centauri 
at the radii of our central HST field. More generally, 
in the central 10 arcmin, the rotation rate is approxi- 
mately linear with radius, V pm t/S pm t ~ R/750". There- 
fore, g w [1 + (-R/750") 2 ] -1 / 2 . These approximations 
are shown as solid lines in the figure. This approxima- 
tion implies that 0.995 < g < 1 for the radii R < 71.7" 
pertaining to the HST sample. 

Given what is already known about proper motion 
rotation from the vdV06 data, our estimates of S pmt 
should be quite accurate. This is despite the fact that 
any solid-body part of the V pm t field cannot be mea- 
sured directly in the HST data. In particular, the value 
of g is insufficient to explain the observed anisotropy 
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Fig. 4. — Observed ratios V pm t/£ pm t (bottom; big dots; ma- 
genta in the on-line version) and g = 2 pmt /2 pm t (top; small dots; 
magenta), as function of projected radius R from the cluster cen- 
ter. The quantities in the nominator and denominator are all av- 
erages over an annulus, and are based on the data of vdV06. The 
solid lines indicate simple approximations to these quantities as 
described in the text. The dashed vertical line indicates the max- 
imum radius for which we use HST proper motions from Paper I. 
The results shown here indicate that our inability to determine 
V pm t from the HST data does not prevent us from accurately de- 
termining the RMS motion S pm t in the tangential proper motion 
direction. 

£pmt/£pmr = 0.983 ± 0.006. Also, the data of vdV06 
show that V\ oa /o~\ os is about a factor two lower than 
V pmt /S pm t. Hence, the influence of rotation is even more 
negligible for the line-of-sight measurements at R = 12" 
of NGB08 (see discussion in Section [5TT]) . 

5.2.4. Gauss-Hermite Moments 

To analyze the shapes of the HST proper motion distri- 
butions we binned the data in concentric circular annuli 
as in Section 15.2.21 For each bin we created separate 
histograms of the proper motions in the pmt and pmr 
directions. We used bins that correspond to 1 kms -1 
for the canonical distance of u Centauri (D — 4.8 kpc; 
vdV06). For each histogram we calculated the best- 
fitting Gaussian and the corresponding Gauss-Hermite 
moments, as in van der Marel et al. (2000). The odd mo- 
ments are zero for distributions that are averaged over 
circles. Figure [5] shows the lowest order non-trivial even 
Gauss-Hermite moments /14 and h 6 . These were calcu- 
lated for annuli that are 8" wide. The average values 
over all radii R < 71.7" are Vpmr = -0.022 ± 0.006, 
Vpmt = -0.024 ± 0.006, h 6tPmi = -0.003 ± 0.006, 
Vpmt = 0.005 ±0.006. The 8th and 10th order moments 
were also calculated, but are not shown here. They were 
found to be consistent with zero, like V For all mo- 
ments, the values for the pmr and pmt directions were 
found to be consistent within the errors. The only mo- 
ment that is non-zero is V with an average value over 
the two proper motion directions of /14 = —0.023 ±0.004. 
This is indicative of proper motion distributions that 
are slightly more flat-topped than a Gaussian, and have 
slightly less extended wings. We will compare this to 
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Fig. 5. — Lowest order even Gauss-Hermite moments /14 (top) 
and ha (bottom) of the observed HST proper motion distributions 
as function of projected radius R from the cluster center. Different 
symbols indicate distributions in the pmr (solid dots) and pmt 
(open circles) directions. Data points in the two directions were 
slightly offset horizontally, for visual clarity. The Gauss-Hermite 
moments are not far from zero, indicating that the proper motion 
distributions are close to Gaussian. There are no systematic trends 
with radius. The curve s sho w predictions of three isotropic models, 
as discussed in Section 16.71 The long-dashed curves (red in the on- 
line version) are for the isotropic core model with no dark mass. 
The solid curves (blue) are for the isotropic cusp model with its 
best-fit IMBH mass A/ B H = 1.8 X 10 4 M Q . Both of these models 
match the observed Gauss-Hermite moments to an average |Afej| < 
0.01. The short-dashed curve (green) is for the isotropic cusp model 
with no dark mass. This model does not fit the observed Gauss- 
Hermite moments. 

model predictions in Section [6.71 

6. DYNAMICAL DATA-MODEL COMPARISON 
6.1. Velocity Anisotropy 

The observed profile of £ pm t/£ pmr shows statistically 
significant deviations from unity, both at small and at 
large radii (see Figure [3b). It is therefore important to 
allow for anisotropy in our models for u> Centauri. 

To illustrate the relation between the model anisotropy 
function /3(r) and the observed quantity £ pm t/£ pmr (i?), 
we show in Figure the predicted E pm t /S pmr (R) for 
models with constant f3 as a function of radius r. Mod- 
els were calculated with {v^^/v 2 ) 1 / 2 ranging from 2/3 
to 3/2, in equal logarithmic steps. At large radii, the 
proper motion anisotropy is a good tracer of the in- 
trinsic anisotropy, consistent with the arguments pre- 
sented in Section l"2. 7. 31 For example, at R = 10 arcmin, 

(Spnt/Epmr) - 1 = T^J^) 1 / 2 - 1], with T = 0.7-0.9 

for the models calculated here 1 and t 0.8 for /3 ra 0. 
At small radii, the behavior of S pm t/S pmr is a more com- 
plicated function of [3. Tangentially anisotropic models 
{(3 < 0) all converge to £ pmt /£ pmr = 1 near the center, 
while radially anisotropic models (0 < (3 < 1) predict 

_ It is evident from Figure [S^ that the observed 
Spmt/S pmr (i?) must tightly constrain the intrinsic (3(r). 



We fitted the combined kinematical data from Figure [3] 
with models that spanned a grid in j3o, /3qo and r a . The 
parameter /3q determines the behavior of £ pmt /E pmr (i?) 
near the center, fioo determines the behavior at large 
radii, and r a is the transition radius. The minimum \ 2 
was obtained for /3 = 0.13 ± 0.02, ^ = -0.53 ± 0.22 
and log r a [arcsec] = 2.86 ± 0.12. The error bars were 
obtained from the x 2 contours of the fit. 

The predicted Ep mt /S pmr (i?) for the best-fit model is 
shown in Figure [Hp. It provides an excellent fit. The 
predictions shown in the figure were calculated for core 
models_with no dark mass. However, the predicted ratio 
£ pmt /£ pmr (i?) depends mostly on /3(r), and is virtually 
indistinguishable when considering instead models with 
a cusp or dark mass. 

The error bars of the fit show that velocity anisotropy 
is detected at high statistical significance in ui Centauri. 
The parameters /3q and Poo imply that the ratio v 2 &n /v 2 
transitions from 0.935 near the center to 1.235 at large 
radii. These are of course not necessarily the values at 
the very center and at infinity, since our model is con- 
strained by proper motion data only over the range 1- 
1000 arcsec. Either way, there is a transition from radial 
anisotropy at small radii to tangential anisotropy at large 
radii. The parameter r a indicates that the transition oc- 
curs at approximately 12 ± 3 arcmin (5 ± 1 core radii). 
This is broadly consistent with the modeling results of 
vdV06. They found (their Section 9.2) that w Centauri 
is slightly radially anisotropic for r < 10 arcmin and that 
it becomes increasingly tangentially anisotropic outside 
this region. 

The centers of globular clusters are generally believed 
to have isotropic velocity distributions because two- 
body relaxation tends to isotropize the orbits. How- 
ever, uj Centauri has a long half-mass relaxation time of 
~ 10 9 96±0 03 years (McLaughlin & van der Marel 2005). 
We showed in Paper I that mass segregation in the cen- 
tral region has not yet progressed to the point of energy 
equipartition. Since full energy equipartition has not yet 
been achieved, it is not surprising that complete velocity 
isotropy has not yet been achieved either. 

A sufficient number of stars is necessary to measure 
the proper motion anisotropy with accuracy. Due to the 
finite number of stars at small radii, the anisotropy is not 
as well constrained at radii R < 10" as it is further out. 
However, because Omega Cen has a very large core, most 
of the stars that are observed near the projected center 
are not actually close to the center in three dimensions 
(see, e.g., Section 6.1 of Paper I). One consequence of this 
is that the projected kinematics predicted near the cen- 
ter are influenced primarily by the anisotropy at larger 
radii. Changing the model anisotropy in the central 10" 
even quite substantially (e.g., ±30% in a r /a t ) does not 
change the model predictions for the HST proper motion 
kinematics by more than a fraction of the error bars, even 
for the innermost data point. Therefore, uncertainties in 
anisotropy near the center do not significantly impact 
our predictions or conclusions. 

6.2. Cusp versus Core Models 

Figures (7^,b show the predicted and observed RMS 
projected velocities for the best-fit (3(r). Each panel 
of this figure shows the data for a\ os (R) and cr pm (i?) 
together in the same panel, with the observed proper 
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Fig. 6. — Dynamical data-model comparisons for the ratio E pm t/E P mr as function of projected radius R from the cluster center. The 
data and the symbols are the same as in Figure[3j:, and are the same in both panels, (a; top panel) Anisotropic models with constant j3 as 

function of radius. The models shown have (v? aIX /v%) 1 ' 2 ranging from 2/3 to 3/2, in equal logarithmic steps. The corresponding j3 values 
are indicated for some of the curves. The horizontal line at 2 pm t/2 pmr = 1 is for the isotropic model. Tangentially anisotropic models 
have higher S pm t /£ p mr than radially anisotropic models. The data are not well fit by models with constant fi. (b) Anisotropic model with 
/3(r) in equation l|15| l chosen the optimize the fit to the data (blue in the on-line version). The model has /3o = 0.13, /3oo = —0.53 and 
r a = 729". The predictions in this figure were calculated for a core model with no dark mass. However, the predicted ratio S pm t/S pmr (iJ) 
depends mostly on j3(r), and is virtually indistinguishable when considering instead models with a cusp or dark mass. 



motions S pm (i?) transformed to km/s using the best-fit 
distance implied by the model. For each model there 
are separate curves predicted for a\ os {R) and a pm (R) 
(dashed and solid curves, respectively), since these quan- 
tities are not the same in an anisotropic model. How- 
ever, the difference is generally smaller than the obser- 
vational error bars over the radial range where both 
types of data are available. There is no line-of-sight 
data available for R < 10", so for visual clarity we 
do not show the predicted a\ os (R) at those radii. At 
small radii ai os (R) exceeds a pin (R) because of the model 
anisotropy, with tri os — cr pm i=a 0.9 km s -1 at R = 10" and 
o'los — 5pm ~ 1-1 kms -1 at R — 1". 

Figure [7^, shows the predictions for the core model fit 
to the projected intensity, and Figure [TJs for the cusp 
model fit. The curves that predict the lowest a at small 
radii (blue in the on-line version) are models with no 
dark mass. The core and cusp models generally predict 
the same dynamics at large radii, but differ slightly at 
small radii. 

It is true in many analytical models with a cusped 
density profile and no central dark mass that the pro- 
jected RMS velocity decreases towards the center (e.g., 
Tremaine et al. 1994). The same can be seen in Figure[7h> 
for our model with 7 = 0.05. Even though the cusp is 
quite shallow in projection, this corresponds to a much 
steeper three-dimensional luminosity density j oc r -1-7 
at asymptotically small radii. Similarly, the decrease in 
RMS velocities towards the center is much stronger in- 



trinsically in three-dimensions than it is in projection. 
Either way, the observed RMS velocities do not show a 
significant dip near the center. In the absence of a dark 
mass, the core models therefore provide a better fit to 
the data than the cusp models. 

The best fits without a dark mass have x 2 = 235.3 and 
249.2 for the core and cusp models, respectively. There 
are 202 data points, yielding JVdf = 197 degrees of free- 
dom (there are five free parameters in the models when 
there is no dark mass) . One expects x ~ ^df ± \ / 2Nuf- 
The anisotropic core model is therefore consistent with 
the data at the 1.9a level, while the cusp model is con- 
sistent at the 2.6cr level. These confidence levels assume 
that all data points in the fit are statistically indepen- 
dent, which is somewhat of a simplification. It is likely 
that samples from different line-of-sight velocity studies 
have some stars in common. However, this cannot be 
fully explored here since Scarpa et al. (2003) and Sol- 
lima et al. (2009) did not list individually which stars 
they observed. 

Much of the excess E{\ 2 ) = x 2 — Nn F > for the 
fits is due to just two data points for a\ os that appear 
anomalously high. The first "discrepant" data point is 
the measurement ai os — 23.0±2.0kms _1 at R = 12 A" by 
NGB08. The second discrepant data point is the average 
<ti os = 19.9 ± 1.0 kms -1 for an annulus centered at R — 
46.0", obtained from the compiled data of vdV06. The 
former contributes 6.8 to the \ 2 of the fit for the core 
model, while the latter contributes 11.3. This adds up 
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Fig. 7. — Dynamical data-model comparisons for the RMS projected velocity as function of projected radius R from the cluster center. 
The data and the symbols are the same as in Figure \3\ and are the same in all panels. For plotting purposes, the proper motions £ pm in 
mas/year were transformed to <r pm in km/s (solid symbols) using equation H25H with the best- fitting distance implied by the models. The 
arrow in the bottom panel indicates how one of the two data points from NGB08 moved due to our improved determination of the cluster 
center in Paper I. Curves are the predictions from various models. Solid curves are for cr pm and short-dashed curves are for O"i os . For visual 
clarity, cti os is not shown for R < 10", where no a\ os data are available. At larger radii, cr pm and rj\ os are quite similar, with differences 
becoming visual only at R > 1000". The curves that predict the lowest a at small radii (blue in the on-line version) are models with no 
dark mass. The remaining curves (red) have some kind of dark mass in the center, (a; top panel) anisotropic models in which /3(r) was 
optimized to fit the data (see Figure [6}, based on the core model for the projected intensity (see Figure |2J. The top curve has the best-fit 
IMBH of mass Mbh = 4.1 x 10 3 Mq. (b; middle panel) Anisotropic models similar to panel (a), but now based on the cusp model for the 
projected intensity (see Figure[2j. The top curve has the best fit IMBH of mass A/bh = 8-7 x 10 3 Mq. (c; bottom panel) Models similar to 
panel (b), but now for an isotropic velocity distribution (which does not fit the data for £ pm t/£ P mr(-f?) in Figure [6}. The quantities cti os 
and <T pm are the same for these isotropic models. The top three curves have different dark masses. The solid curve is for the best-fitting 
IMBH mass M BH = 1.8 X 10 4 Mq. The dotted curve is for the IMBH mass M B H = 4.0 X 10 4 Mq advocated by NGB08. The long-dashed 
curve is for the best fitting dark cluster, which has M^ark = 2.0 X 10 4 Mq and a(j ar k = 3.0". It follows from analysis of the predictions 
in this figure, and in particular the top panel, that the presence of an IMBH in u) Centauri is not required, although IMBHs with masses 
A/bh i; 1.2 X 10 4 Mq cannot be ruled out at ~ 1<t confidence. 



to 18.1, which is almost half of E(\ 2 ) = 38.3. 

The discrepant NGB08 data point pertains to similar 
radius as the other field that they observed, which is at 
R = 12.1" from the cluster center (see Paper I). The 
measurement <ji os = 18. 6± 1.6kms _1 for that latter field 
is more consistent with the other measurements in our 
combined data set from Sectional The two NGB08 mea- 
surements are inconsistent with each other at the 1.7 a 
level, despite the fact that that they pertain to the same 
radius, and despite the fact that the proper motions in 
these fields are consistent with each other (Paper I). We 



believe that this may be due to underestimates in the 
NGB08 error bars, since they did not explicitly include 
the contribution of shot noise (from the finite number 
of stars) to their error bars. Similarly, the measurement 
in the vdV06 annulus at R = 46.0" appears inconsistent 
with other measurements in that same data set. The 
weighted average of the two adjacent annuli is lower and 
inconsistent with it at the 2.1a level. This suggests that 
both discrepant measurements that contribute dispro- 
portionately to x 2 may be spuriously high for unknown 
reasons. 
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Fig. 8. — The \ 2 of the model fits to the kinematic data in 
Figure[3] as function of IMBH mass Mbh- The resulting parabolae 
are, from bottom left to top right, respectively, obtained for the 
following models: anisotropic models with a core (black in the on- 
line version; see Figure[7^ for actual model predictions) ; anisotropic 
models with a cusp (red; see Figure[7p); and isotropic models with 
a cusp (blue; see Figure [7};). Centered above the minimum of each 
parabola is a triangle with an associated error bar that indicates the 
best fit IMBH mass and its formal 1<t uncertainty. The anisotropic 
core models provide the overall best fits to the kinematics. The 
IMBH that optimizes the \ 2 f° r those models is not statistically 
significant. The small difference in \ 2 (vertical height) between the 
parabolae for the anisotropic core and cusp models is not a reason 
to prefer to core models over cusp models, in particular because 
the latter actually provide a somewhat better fit to the photometry 
(see discussion in the text). 

Evidently, there are some internal discrepancies in the 
datasets themselves. Hence, the observed values of E(x 2 ) 
should not be taken as a sign of shortcomings in the mod- 
els. In fact, the values of E{x 2 ) are remarkably low, given 
that we are fitting kinematical data compiled from 10 
different original sources. Potential uncorrected system- 
atics between different studies could easily have induced 
much larger data-model discrepancies. 

One clear result from Figures [7K,b is that the data 
are well fit by the models out to the largest available 
radii. This contradicts arguments put forward by Scarpa 
et al. (2003) that the large radii kinematics of u> Cen- 
tauri may be inconsistent with traditional theories of 
gravity. Our conclusion on this issue agrees with that of 
McLaughlin & Meylan (2003) and Sollima et al. (2009). 

6.3. IMBH 

To test for the presence of an IMBH we have also con- 
structed models with Mbh > (and the same f3(r) pro- 
file). Models with no IMBH have one free parameter 
less than models with an IMBH. One would therefore on 
average expect them to fit worse by A\ 2 = 1- More- 
over, random variations can add another Ax 2 = 1 (at 
la confidence), even when the no-IMBH model is cor- 
rect. Therefore, an IMBH is statistically significant in 
the model only if it improves the fit over a no-IMBH 
model by A% 2 > 2. 

Figure [8] shows the curves of x 2 (Mbh) that quantify 



the quality of the fit. For the core model, the best-fit 
IMBH has mass Af B H = 4.1 x 10 3 M Q . This improves 
the fit over the no-IMBH model by only A% 2 = 1.8. This 
is not statistically significant. For the core model we can 
therefore determine only an upper limit to the mass of 
a possible IMBH, namely M B h < 7.4 x 10 3 M Q at Ict 
confidence. Figure [7^, compares the best-fit core models 
with and without an IMBH. The inclusion of the IMBH 
makes very little difference to the predictions. The dif- 
ference for the innermost data point at R = 1.1" is only 
0.8 kins" 1 , which is less than a quarter of the observa- 
tional error bar at that radius. 

For the cusp model, the inclusion of an IMBH im- 
proves the fit by A% 2 = 9.8, which is significant at 
the 3. Ict level. The best-fit IMBH mass is M B h = 
(8.7±2.9) x 10 3 M Q (with the error bar determined by the 
requirement Ax 2 < 1). Figure [TJa compares the best-fit 
cusp models with and without an IMBH. The inclusion 
of the IMBH provides a small but visible improvement 
in the fit near the center. The reason that the IMBH is 
more massive and significant in the cusp model is due to 
the dip in RMS velocity predicted by these models when 
there is no dark mass. 

In assessing the meaning of formal error bars on Mbh, 
it is important to keep in mind what one may call "black 
hole bias". Statistical error bar estimates take into ac- 
count only the random Gaussian uncertainties in the 
data. They do not take into account the systematic resid- 
uals that are inevitably present in any astronomical data- 
model comparison. In the present context these may in- 
clude oversimplifications in the model assumptions (e.g., 
sphericity, equilibrium, the adopted density parameteri- 
zation, no mass segregation among the luminous stars), 
the data reduction process (e.g., inability to measure ab- 
solute proper motion rotation), or the meaning of the 
data themselves (e.g., inability to separate orbital mo- 
tion from potential binary motion). By adding Mbh 
as an extra free parameter to the fit, the model gains 
in ability to fit these systematic residuals that are not 
actually due to an IMBH. Because masses are positive 
definite (A/bh > 0), this generally tends to bias in the 
direction of invoking more mass then there actually is, 
as compared to the opposite. Hence, Mbh tends to be 
biased high. Marginally significant non-zero Mbh values 
with no tell-tale signs of an IMBH (e.g., fast-moving stars 
or an i? -1 / 2 increase in a) may simply be indications of 
small systematic issues in the data-model comparison un- 
related to an IMBH. 

The core model with no IMBH fits the kinematical 
data a little better than the cusp model with an IMBH 
(see Figure [5]). The difference is A\ 2 = 4.0, which is 
nominally significant at the 2. Oct level. However, as dis- 
cussed in Section [U the cusp model nominally provides 
a better fit to the photometric data. The difference for 
a generalized nukcr fit is A\ 2 = 6.8, which is nominally 
significant at the 2.6ct level. One could in principle con- 
struct a grand-total A% 2 that sums the photometric and 
kinematic values. However, we have not done this here. 
The arguments about whether a core or cusp model bet- 
ter fits the photometry are complex (see Section[4|) . They 
are not easily captured in a single number without hav- 
ing to make arbitrary choices about, e.g., what radii to 
include in the comparison. Either way, when comparing 
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both the photometric and kinematic information for the 
core and cusp models we feel that there is no sound sta- 
tistical basis to prefer one model over the other. In view 
of this, we conclude that the existing data do not imply 
that an IMBH is present in uj Ccntauri. However, IMBHs 
with masses Mbh ^ 1-2 x 10 4 Mq cannot be ruled out 
at la confidence (or < 1.8 x 10 4 Mq at 3cr confidence). 

6.4. Distance and Mass-to-Light Ratio 

The best-fit distance and mass-to-light ratio depend 
only slightly on the particular model we adopt. The core 
model without a dark mass has D — 4.70 ± 0.06 kpc 
and Ty = 2.64 ± 0.03 (here and henceforth, mass-to- 
light ratios are given in V-band solar units). The cusp 
model with an IMBH has D = 4.75 ± 0.06 kpc and Ty = 
2.59 ±0.03. 

Our results for D and Ty are entirely consistent with 
those inferred by vdV06 from fitting only ground-based 
projected intensity and kinematical data. They found 
that D = 4.8 ± 0.3 kpc and Ty = 2.5 ± 0.1. Their 
modeling technique was considerably more sophisticated 
than ours, allowing for axisymmetry, inclination, clliptic- 
ity variations with radius, and arbitrary profiles of rota- 
tion and velocity-dispersion anisotropy. The reason that 
our error bars are smaller than theirs is likely a combi- 
nation of two effects. First, we have considerably more 
data at our disposal than did vdV06, and in particular we 
have a factor ~ 11 more stars with measured proper mo- 
tions. This produces a significant decrease in the random 
uncertainties. But second, in our simpler models we can- 
not explore the full variety of phase-space distributions 
that might fit the data with non-standard distances or 
mass-to-light ratios. Therefore, we are not quantifying 
systematic errors as rigorously as vdV06, and as a result 
our errors are artificially lowered compared to theirs. 

The fact that our D and Ty agree with those of 
vdV06 indicates that there is no reason to mistrust that 
our spherical models can address the central mass dis- 
tribution of uj Ccntauri with credible accuracy. This 
is further supported by the facts that both the elon- 
gation (Geyer, Nelles, & Hopp 1983) and the rotation 
rate (Section I5.2.3|) decrease towards the cluster center. 
To get good results we did have to make sure to model 
RMS velocities (which sum rotation velocities and veloc- 
ity dispersions in quadrature) and not merely velocity 
dispersions. Also, one should fit only quantities that are 
properly averaged over annuli on the projected plane of 
the sky. If these things are not done, then biased esti- 
mates are obtained for D and Ty, as demonstrated in 
Appendix C of vdV06. 

6.5. Isotropic Models 

Isotropic models predict that E pm t/S pmr = 1 at all 
radii, so they are not generally appropriate for uj Cen- 
tauri. However, isotropic models have the advantage that 
the full velocity distributions £i so (c, R) can be calculated 
as described in Section 12.51 Also, isotropic models are 
often used for globular clusters, so it is useful to under- 
stand how the predictions of such models deviate from 
more general anisotropic models. 

Figure[8]shows the curve of x 2 (Mbh) for isotropic mod- 
els with an IMBH. The models fit significantly worse 
than the anisotropic models, as shown by the higher 
X 2 . Within the realm of isotropic models, the best fit 



is provided by M B h = (1-8 ± 0.3) x 10 4 M . This model 
has D = 4.81 ± 0.06 kpc and Ty = 2.61 ± 0.03. Fig- 
ure [7]; shows the data-model comparison for the isotropic 
model with this IMBH. We also show the predictions for 
Mbh = 0, as well as for the value Af B H = 4.0 x 10 4 M Q 
advocated by NGB08. 

The NGB08 IMBH mass was based on their measure- 
ment of CTi os = 23.0 ± 2.0 kms -1 at a position they be- 
lieved to be the cluster center. This high IMBH mass 
is now clearly ruled out, even for an isotropic model, 
for two reasons. First, Paper I showed that this mea- 
surement actually applies to a larger distance R = 12.4" 
from the cluster center. The arrow in Figure [7b demon- 
strates how this moved the corresponding line-of-sight 
velocity measurement. And second, the high measured 
RMS velocity was not confirmed by proper motion mea- 
surements, either at the same position or on the actual 
cluster center. Note in Figure^ that the isotropic model 
with the NGB08 IMBH does appear to match well the 
data at R = 5-7". However, it overpredicts the data at 
both R < 4" and R = 8-10". Since all data points were 
derived in the same manner and are equally valid, there 
is no reason attribute more weight to the data at R = 5- 
7". In particular, the data at R < 4" contain most of 
the information on the very central mass distribution, 
and they clearly rule out an IMBH mass as high as ad- 
vocated by NGB08 (independent of the model anisotropy 
near the center). 

It is evident from Figure [8] that isotropic models re- 
quire a higher IMBH mass than the best-fit anisotropic 
models. This is not generally true for isotropic mod- 
els, but is a specific consequence of the fact that uj Cen- 
tauri is radially anisotropic near the center and tangen- 
tially anisotropic at large radii. This affects the general 
gradient in RMS velocity from small to large radii (see 
e.g., Figures 10 and 11 of van der Marel 1994). The 
anisotropic model without a dark mass has a steeper gra- 
dient towards the center than the corresponding isotropic 
model. The anisotropic model therefore requires less 
dark mass to fit the observed gradient. This indicates 
that one should be careful in application of isotropic 
models to studies of IMBHs. Depending on the exact 
anisotropy of the cluster, the assumption of isotropy can 
lead to either an overestimated or an underestimated 
IMBH mass, or it can require an IMBH when none is 
present. 

6.6. Dark cluster 

When a dark mass is clearly required by the data, as 
in the case of isotropic models with a cusp, one can ask 
whether an extended dark cluster may provide a better 
fit than an IMBH. Figure [9] shows a two-dimensional con- 
tour plot of A% 2 in the parameter space of (adark, -Miark) 
for the isotropic cusp model. The best fit is obtained for 
-Mdaric = 2.0 x 10 4 Mq and Odark = 3". The predictions 
of this model are shown in Figure [7J:. It has a somewhat 
shallower increase in RMS velocity towards the center 
than the IMBH model. The dark-cluster model makes an 
improvement of only A% 2 = 1.7 over the IMBH model, 
while increasing the number of free parameters by one. 
This is not statistically significant. Therefore, one can 
interpret Figure [9] to mean that if there is a dark mass 
in uj Centauri, then its extent must be < 7" at la confi- 
dence. 
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Fig. 9. — Contours of A% 2 in the two-dimensional plane of 
( a darki A^dark) f° r the fit of isotropic cusp models to the data in 
Figure [3] The distance D and mass-to-light ratio Ty were op- 
timized separately at each (adark i -Mlark) f° r which a model was 
calculated (indicated with small dots). The model for which the 
minimum \ 2 was obtained is indicated by the big dot. The first 
three contours (solid) around the minimum correspond to la, 2a, 
and 3cr (bold) confidence contours, respectively. Each subsequent 
contour (dashed) increases A\ 2 by a constant factor. Models with 
a central dark points mass (an IMBH; i.e., adark = 0) are shown 
on the vertical axis on the left. Models with no central dark mass 
lie on the horizontal axis on the bottom. If there is a dark mass 
in ui Centauri (as must be the case in isotropic models, but these 
do not fit the data for 2 pm t/2p mr (-R) in Figure[6]l then its extent 
must be < 7" at la confidence. 

The central density of the best-fitting dark cluster is 
/Oo,dark = Af dark /(47ra^ ark /3), according to equation 
Using equation ^ to transform angular units to physi- 
cal units, this yields for the core model po.dark = 1-4 X 
10 7 Mq pc -3 . The central density of luminous mat- 
ter in this model is only po,ium = 3.3 x 10 3 Mq pc~ 3 . 
This implies a central M/L = Ty[l + (po,dark//Oo,ium)] 
1 . 1 x 10 4 . This applies to an isotropic model with 
-^dark = 2.0 x 10 4 Mq. For an anisotropic model with a 
lower A/dark, the required M/L scales approximately lin- 
early with A/dark- The (maximum) size adark of the dark 
mass fit is relatively insensitive to the exact anisotropy. 
Either way, if u> Centauri does indeed have a dark cen- 
tral cluster, this would be indicative of quite an extreme 
amount of mass segregation or a very top-heavy initial 
mass function. If a dark mass is invoked to explain the 
data, then an IMBH may be easier to explain than such 
a concentrated dark cluster. 

6.7. Velocity Distribution Shapes 

For our isotropic models we calculated the projected 
velocity distributions C[ SO (v,R) as described in Sec- 
tion I2.5[ an d from this the Gauss-Hermite moments 
hi.iso(R)- We did not take the effect of observational 
proper motion uncertainties into account in the models 
(essentially a convolution of the model predictions with 
the error distribution). This has negligible effect on the 
Gauss-Hermite moments because the observational er- 
rors are much less than the intrinsic cluster dispersion. 
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Fig. 10. — Symmetrized histograms of the observed proper 
motion velocities for apertures of radii R = 3" (bold; red in the 
on-line version) and R = 10" (thin; blue), respectively. Proper 
motions in both the pmr and pmt directions are included. They 
were transformed to physical units of km/s using the best-fit model 
distance. Typical error bars are shown on the left; they were offset 
horizontally from the histograms for visual clarity. The predictions 
for the isotropic cusp model with the best fit IMBH mass Mbh = 
1.8 X 10 4 Mq are over-plotted as curves (solid for the R = 3" 
aperture, and dashed for the R = 10" aperture). The inset on 
the right shows a blow-up of the profile wings on an expanded 
vertical scale. The models predict broader wings than observed, 
but the difference is not statistically significant (as discussed in 
the text). This is because the total area under the model curves 
beyond 60kms — 1 corresponds to less than 1 predicted star. 



The curves over-plotted in Figure [5] show the model pre- 
dictions for three models: the isotropic core model with 
no dark mass, the isotropic cusp model with its best-fit 
Mbh = 1.8 x 10 4 Mq, and the isotropic cusp model with 
no dark mass. 

The core model with no dark mass and the cusp model 
with an IMBH both fit the observed Gauss-Hermite mo- 
ments within the error bars. Over the range 5" < R < 
71.7", the predicted Gauss-Hermite moments for these 
models are approximately constant. For the core model 
without dark mass they are h4^ so (R) « —0.013 and 
he,iso(R) ~ —0.002, and for the cusp model with an 
IMBH they are h^ iso (R) » -0.034 and h 6 ,i so (R) ~ 0.009. 
These averages do deviate in statistically significant man- 
ner from the observed averages, hi = — 0.023± 0.004 and 
h 6 = 0.001±0.004 (see Section I5T4"]) . This is not surpris- 
ing, given that we have already established that uj Cen- 
tauri does not have an isotropic velocity distribution, not 
even near its center. However, the differences between 
the predicted and the average observed Gauss-Hermite 
moments are quite small \Ahi\ < 0.01. For comparison, 
models that range in anisotropy from fully tangential to 
fully radial can easily span the range Ahi ss ±0.2 (e.g., 
van der Marel & Franx 1993). The measured Gauss- 
Hermite moments are therefore consistent with the fact 
that the velocity distribution is not far from isotropic 
near the center of lu Centauri. 

The cusp model without dark mass does not fit the 
observed Gauss-Hermite moments. As discussed in Sec- 
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tion 16.21 this model has a strong decrease in its intrinsic 
(unprojected) RMS velocity towards the center. The de- 
crease in the projected RMS velocity is relatively mild, 
as shown in Figure [7J:. However, there is a strong gra- 
dient in the RMS velocity along the line of sight. This 
leads to strongly non-Gaussian velocity distributions. In 
a central aperture of R — 8" radius, as used for the 
innermost observational data point in Figure [3 the pre- 
dictions are h^o — 0.063 and 7i6,iso = —0.063. This 
deviates for each order by ~ 3<r from the observed val- 
ues h 4 = -0.026 ± 0.029 and h 6 = 0.031 ± 0.029 (aver- 
aged over the pmr and pmt directions). While this result 
was derived for isotropic models, the large predicted de- 
viations from Gaussian velocity distributions are likely 
generic, and would exist in anisotropic models as well. 
So if lo Centauri has a central cusp in its projected in- 
tensity profile, then the observed velocity distributions 
rule out models without a dark mass. Of course, Fig- 
ure [7b, c show that such models also do not provide good 
fits to the projected profiles of RMS velocity. 

Figure [5J shows that an isotropic model with an IMBH 
predicts increased /14 and he in the central R < 5", cor- 
responding to profiles with broader wings than a Gaus- 
sian. To do a more detailed data-model comparison for 
the wings we show in Figure [10] the histograms of the 
observed proper motions for apertures of radii R = 3" 
and R = 10". These apertures contain N = 43 and 
N = 585 stars with well-measured proper motions, re- 
spectively. The proper motion coordinates in the pmr 
and pmt direction were both included in each histogram, 
as appropriate for an isotropic velocity distribution. The 
histograms were symmetrized to decrease the uncertain- 
ties. The predictions for the isotropic cusp model with 
M BH = 1-8 x 10 4 Mq are over-plotted as curves. The 
predicted velocity distribution for the R = 3" aperture 
has more extended wings than the observed histogram. 
The wings are much reduced in the predicted velocity 
distribution for the R = 10" aperture. This is because 
none of the added stars at 3" < R < 10" are close to 
the black hole, so they do not move at unusually high 
velocities. 

The observed histograms do not extend beyond 
60kms . The maximum observed value of either |u pm t| 
or |u pmr | inside R = 10" is 59.4kms _1 . No stars were re- 
jected as non cluster members inside this radius because 
of their fast motion (see Section [5.2.1[) . To assess whether 
these observations are statistically consistent with the 
wings of the predicted distributions, let / be the fraction 
of the normalized model distribution at |«| > 60kms . 
According to the models, this is / = 3.8 x 10~ 3 for the 
R = 3" aperture and / = 3.9 x 10~ 4 for the R = 10" aper- 
ture. The expectation value for the number of stars with 
a velocity component at \v\ > 60 kms^ 1 is E = 2Nf. 
The probability that all of the 27V measured velocities 
have \v\ < eOkms" 1 is P = (1 - f) 2N . This yields for 
the R = 3" aperture that E = 0.32 and P = 0.72, while 
for the R = 10" aperture E = 0.46 and P = 0.63. 

Since less than a single fast-moving star was predicted, 
the fact none were observed is not statistically inconsis- 
tent with the models. This does not change when we 
use larger apertures, since the models do not predict 
any fast-moving stars outside R = 10". It also does not 
change by using smaller apertures, since in that case the 



sample of observed stars becomes very small. Therefore, 
the absence of fast-moving stars in the data cannot be 
used to independently rule out the presence of an IMBH 
with mass Mbh = 1.8 x 10 4 M©. To increase the pre- 
dicted number of fast-moving stars one would need to 
increase the number of stars with measured proper mo- 
tions, i.e., extend the kinematical observations to fainter 
magnitudes. We have experimented with looser cuts on 
our proper motion catalog. This allows us to increase 
the sample size by a factor ~ 2 (see Paper I). In this 
extended catalog there is no star with a proper motion 
coordinate in excess of 65 kms -1 within R = 10". This 
provides somewhat tighter constraints on the presence of 
an IMBH, with the quoted values of P roughly decreasing 
to their squared values. Even then, it is still not possible 
to rule out an IMBH of mass M B h = 1.8 x 10 4 M© with 
reasonable confidence. 

Alternatively, E increases and P decreases when the 
BH mass is increased. For the IMBH of mass Mbh = 
4.0 x 10 4 Mq suggested by NGB08, and sticking as be- 
fore with our high-quality proper motion catalog of Sec- 
tion [5XT1 we obtain E = 2.9 and P = 0.05 for the 
aperture with R = 10". This implies that the ab- 
sence of fast-moving stars rules out all IMBHs with 
M B h > 4.0 x 10 4 M© at 95% confidence. This result is not 
very competitive for lo Centauri, since we already know 
from modeling of the RMS velocities that any IMBH 
must have a mass at least 3.3 times lower than this any- 
way. However, this method has the advantage that it 
does not require any accurate measurement of RMS ve- 
locities. It may therefore be useful for other clusters, in 
cases when the observational errors may be too large for 
accurate measurement of the RMS velocity profile. 

6.8. Brownian Motion 

We have assumed in the modeling that any IMBH must 
reside at the cluster center. In reality, an IMBH would 
exhibit Brownian motion as a result of its energy equipar- 
tition with surrounding stars. Chatterjee, Hernquist & 
Loeb (2002) derived for the one-dimensional RMS offset 
from the cluster center that 

xrms = (2m/9M BH ) 1/2 r . (34) 

This was derived for an analytical Plummer model with 
core radius ro, which provides a reasonable descrip- 
tion of (non core-collapsed) globular clusters. Chat- 
terjee et al. considered single-mass models with stars 
of mass fh. However, their results should be valid for 
multi-mass models as well, provided that the quantity 
fh = (m 2 )/(m) is defined appropriately in terms of av- 
erages over the stellar mass function (Merritt, Berczik & 
Laun 2007). 

We take for ro the King core radius of 141" (McLaugh- 
lin & van der Marel 2005) and for Mbh the upper limit 
of 1.2 x 10 4 Mq derived in Section [631 For a Salpeter 
mass function from 0.1-10 M©, fh = 1.27 M©. We then 
obtain xrms = 0.69". This is significantly less than the 
1.55" radius of the smallest aperture used in the kinemat- 
ical analysis (see Section 15.2.21) . Therefore, the expected 
Brownian motion of an IMBH should not affect signifi- 
cantly the analysis that we have presented here. 

IMBHs with larger masses than those considered here, 
placed at considerable offsets from the cluster center 
(consistent with their expected Brownian motion), might 
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TABLE 1 

Globular clusters with IMBH constraints 



cluster 


D 


Vb 


M/L v 


Aftot 


o" e 


M BH [Tre] 


Af B H [obs] 


data 


A/ BH [obs] / Af tot 




(kpc) 


(mag) 


(©) 


(M e ) 


(km/s) 


(M e ) 








(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


47 Tuc 


4.5 6 


3.91 6 


2.4 6 


1.1 x 10 6 


9.8 fc 


7.3 x 10 2 


< 1.5 x 10 3s 


discrete PM 


< 0.13% 


Omega Cen 


4.7° 


3.14 6 


2.6 a 


2.8 x 10 6 


15.7 J 


4.8 x 10 3 


< 1.2 x 10 4a 


discrete PM+LOS 


< 0.43% 


M15 


10.3' 


5.89 e 


1.6 C 


6.2 x 10 5 


12.1 c 


1.7 x 10 3 


< 4.4 x 10 3h 


discrete PM+LOS 


< 0.71% 


NGC 2298 


10.7 6 


9.36 6 


1.9 6 


3.4 x 10 4 


2.4 fe 


2.6 x 10° 


< 3.4 x 10 2 f 


mass segregation 


< 1.0% 


Gl 


770 1 


13.49* 


2.8 d 


5.7 x 10 6 


25.1* 


3.2 x 10 4 


(1.8 ± 0.5) x 10 4d 


unresolved LOS 


(0.32 ± 0.09)% 



note. — Column (1) lists the cluster name. Column (2) lists the cluster distance D. Column (3) lists the extinction-corrected 
total apparent magnitude Vo- Column (4) lists the V-band mass-to- light ratio M/Lv in units of Mq/ Lq v . Column (5) lists 
the total cluster mass M to t derived from columns (2)-(4). Column (6) lists the effective velocity dispersion a c . Column (7) 
lists the value of Mbh predicted by the A/bh vs. cr e relationship for galaxies of Tremaine et al. (2002), using the value of a e 
in column (6). The uncertainty in the predicted Mbh due to the uncertainties in the relation itself is a factor ~ 1.5. There 
is also intrinsic scatter about the relationship of a factor ~ 2. Column (8) lists the observational constraint on the mass 
Mbh of any IMBH. Column (9) lists the type of data on which the observational constraint in column (8) is based. This 
can be either discrete velocities in the proper motion (PM) and/or line-of-sight (LOS) directions; unresolved observations of 
line-of-sight kinematics; or measurements of the amount of mass segregation, which can be compared to model predictions. 
Column (10) lists the observational constraint on the ratio MBH/M to t- Sources of the numbers in the table are: a This paper. 
6 McLaughlin & van der Marel (2005). c Gerssen et al. (2002). d Gebhardt et al. (2005). e Harris (1996). 1 Pasquato et 
al. (2009). 9 McLaughlin et al. (2006). h Gerssen et al. (2003). * Meylan et al. (2001). 3 Line-of-sight velocity dispersion in 
an aperture of size equal to the effective radius, calculated from the models in this paper; this is 85% of the central value. 



The central model velocity dispersion value in McLaughlin 
derived in footnote j. 1 van den Bosch et al. (2006). 



van der Marel (2005), multiplied by the same 85% correction 



have gone unnoticed in our analysis. However, there is 
no independent motivation to consider such a situation. 
In Paper I we showed that the kinematical center of stars 
in to Centauri agrees with the star count center to within 
its ~ 2" uncertainty. 

7. GLOBULAR CLUSTER IMBH DEMOGRAPHICS 

The self-gravity of a cluster gives it a negative heat 
capacity that makes it vulnerable to the so-called 
"gravothermal catastrophe" . This makes the core col- 
lapse on a timescale proportional to the half-mass relax- 
ation time. This is normally halted by formation and/or 
heating of binaries. However, the most massive stars in 
a cluster tend to undergo core collapse more or less inde- 
pendently of the other cluster stars on a timescale that 
is less than the core collapse time for the cluster as a 
whole. If this happens faster than the main sequence 
evolution time of these stars (i.e., before they become 
compact remnants) , then they will undergo physical col- 
lisions that lead to the formation of a single Very Massive 
Star (VMS). Such a VMS could plausibly evolve into an 
IMBH. 

Portegies Zwart & McMillan (2002) proposed the 
aforementioned scenario and performed iV-body simula- 
tions to show that a VMS might grow to MvMs/-^tot ~ 
0.1%, with M to t being the total cluster mass. Gurkan 
et al. (2004) and Freitag et al. (2006) used Monte-Carlo 
simulations with prescriptions and models for stellar col- 
lisions and found that the resulting VMS could even grow 
up to MvMs/-Mtot ~ 1%- However, it is not clear what 
the resulting IMBH mass might be. The evolution of a 
VMS may yield significant mass loss before the final col- 
lapse. Also, a VMS will only form through this scenario 
in clusters with very short initial half-mass relaxation 
times, < 30 Myr. Some star clusters (e.g., the Arches 
cluster) are known to have such short relaxation times. 



However, most globular clusters have much longer cur- 
rent relaxation times, although it is not straightforward 
to estimate what their relaxation time at birth may have 
been (Ardi et al. 2008). 

Alternative scenarios for IMBH formation in globular 
clusters have been proposed as well (see e.g. the review in 
van der Marel 2004) . This includes the merging of stellar- 
mass black holes. O'Leary et al. (2006) found that this 
may lead to growth up to MiMBH/M to t ~ O-l^o- Whether 
this indeed happens depends on the typical BH merger 
recoil speed. Significant recoils will kick black holes out 
of the cluster before they have a chance to grow. 

Overall, the predictions for IMBH formation in glob- 
ular clusters have significant uncertainties. Observa- 
tions will therefore need to be the primary guide for the 
study of IMBH demographics in globular clusters. For 
to Centauri one may combine the distance and mass-to- 
light ratio derived in Section 16.41 with the extinction- 
corrected total V-band magnitude to determine that 
M tot = 2.8 x 10 6 M . The IMBH mass upper limit from 
Section |6~31 therefore corresponds to M-Q-a/Mtot ^ 0.43%. 
Table [T] lists the corresponding constraints for the other 
globular clusters that have been well-studied dynami- 
cally, as discussed in Section [1] We also include the limit 
recently reported by Pasquato et al. (2009) based on a 
study of the observed mass segregation in the cluster 
NGC 2298. An IMBH is expected to reduce the amount 
of mass-segregation that is otherwise expected in a well- 
relaxed cluster. The suggested IMBH in Gl weighs in at 
-Mbh /Mtot = (0.32 ± 0.09)%. The four other clusters in 
the table only provide upper limits in the range 0.1-1%. 
Therefore, the data allow us to probe in a mass range 
that is of interest in view of existing theories. 

The IMBHs suspected in some globular clusters may 
have been similar to the initial seeds that grew to super- 
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massive black holes (SMBHs) in the centers of galaxies. 
It is therefore of interest to compare the results in Ta- 
ble Q] to what is known about SMBHs. The masses of 
SMBHs scale with either the galaxy bulge mass Mb u i, 
or alternatively, with the fourth power of the veloc- 
ity dispersion a. Haering & Rix (2004) found that 
M BH /M hvL i w 0.14 ± 0.04%. This is below where it has 
been possible to probe for most globular clusters. Only 
for 47 Tuc is there a limit M Bii /M tot < 0.13% (McLaugh- 
lin et al. 2006). Tremaine et al. (2002) discussed the 
Mbh~(J relation proposed for SMBHs. In Tabic [T] wc 
list for each cluster the value of Mbh predicted by their 
relation. All of the available Mbh upper limits are at 
least a factor of 2 above these predictions. Therefore, the 
data are insufficient to confirm or rule out the hypothesis 
that globular clusters generally have IMBHs that follow 
the same relations as established for SMBHs. However, 
for the case of Gl the suggested IMBH mass does seem 
reasonably consistent with these relations (Gebhardt et 
al. 2005). 

8. DISCUSSION AND CONCLUSIONS 

We have presented a detailed dynamical analysis of 
the star-count, surface-brightness, and kinematical data 
available for u> Centauri, with a particular focus on the 
new HST data presented in Paper I. Based on the ob- 
served profile of the projected density, our models use the 
Jeans equation to yield predictions for the projected pro- 
files of RMS linc-of-sight velocity a or proper motion £ 
as function of projected distance R from the cluster cen- 
ter, in each of the three orthogonal coordinate directions 
(line of sight, proper motion radial, and proper motion 
tangential). We do not include the effect of mass segre- 
gation on the model predictions at large radii, but we do 
take into account the observed (partial) energy equipar- 
tition between stars of different masses. The predictions 
are compared in a x 2 sense to observations, to infer the 
best-fit values of the model parameters and their error 
bars. 

The model parameters arc uniquely determined by 
the data. The profile (3(r) of the intrinsic velocity 
anisotropy is tightly constrained by the observed profile 
[S pm t/S pmr ](i?). The models are therefore unaffected by 
the mass-anisotropy degeneracy that often plagues line- 
of-sight velocity studies. The mass-to-light ratio T of 
the stellar population is fixed by normalization of the 
observed a\ os (R) and £ pm (i?) at large radii. The dis- 
tance D is set by comparison of the observed £ pm (i?) 
in mas/yr to the predicted <7 pm (r) in km/s. And the 
presence and properties of any dark mass in the center, 
such as an IMBH or dark cluster, are constrained by the 
behaviors of a\ os {R) and £ pm (i?) at small radii. 

The new HST proper motion data provide a significant 
improvement in our ability to constrain the dynamical 
models of lu Centauri, as compared to previous model- 
ing efforts (e.g., vdV06; NGB08). The HST sample in- 
creases the number of stars with one or more accurately 
measured velocity components by about an order of mag- 
nitude. This allows for the first time the measurement 
of velocity dispersions at R < 10". We have derived 
kinematical profiles down to R ~ 1", at which point 
there cease to be enough stars with measured proper 
motions to accurately measure a dispersion. The pro- 
file of E pm (i?) is consistent with being flat in the central 



R < 15". The absence of an increase in RMS stellar ve- 
locities towards the center provides strong constraints on 
the possible presence of any dark mass. 

The HST star count data from Paper I provide an up- 
per limit 7 < 0.07 to the central logarithmic cusp of 
slope of the projected density profile. The innermost 
data points are well matched by previously published 
King and Wilson models with a flat core. A generalized 
nuker profile fits best at all cluster radii when a shallow 
cusp is invoked, but this overpredicts the innermost dat- 
apoints at R < 4". Dynamical model predictions depend 
sensitively on whether or not there is a density cusp, so 
in our analysis we considered separately models with a 
core (7 = 0) and models with a cusp (7 = 0.05). 

The inferred intrinsic anisotropy profile, distance, and 
stellar mass-to-light ratio are all insensitive to the exact 
light and mass distributions near the center. We find 
that the anisotropy changes from slightly radial near the 
center to significantly tangential in the outer parts. The 
presence of velocity anisotropy is not surprising, given 
the long two-body relaxation time of u Centauri and the 
existence of other evidence that the cluster is not yet 
fully relaxed. The inferred distance is 4.73 ± 0.09 kpc 
and the mass-to-light ratio is Yy = 2.62 ± 0.06. The 
error bars take into account both the uncertainties in the 
kinematical data and the uncertainty in the cusp slope 7. 
These results match extremely well with those obtained 
by vdV06 using more sophisticated modeling tools. This 
agreement supports our assertion that spherical models 
based on the Jeans equation are sufficient for the goals 
of the present paper, and in particular for constraining 
the mass distribution near the center. The facts that the 
ellipticity and rotation of to Centauri are both known to 
decrease towards the center provide further motivation 
for this assertion. 

Models with a core provide a good fit to the kine- 
matical data without any dark mass. Since a core is 
also consistent with the observed density profile, this im- 
plies that the presence of an IMBH is not required in 
u> Centauri. By contrast, models with a shallow cusp 
provide a good fit only when an IMBH is invoked, with 
A/bh = (8.7 ± 2.9) x 10 3 M Q . Indeed, an IMBH can in- 
duce a shallow density cusp (Baumgardt, Makino & Hut 
2005). Without a dark mass, the predicted RMS veloci- 
ties in the cusp model decrease towards the cluster center 
and the predicted velocity distributions become highly 
non-Gaussian. Neither prediction is consistent with the 
observations. 

Upon combination of the predictions for the core and 
cusp models we obtain as final result an upper limit to 
the mass of any possible IMBH of 1.2 x 10 4 M© at la 
confidence (or 1.8 x 10 4 M© at 3a confidence). The la 
limit corresponds to MBn/M to t <j 0.43%. This is of sim- 
ilar magnitude as the range that has now been probed 
for several other clusters as well. However, lower values 
of MBu/Mtot will need to be probed to definitively assess 
the most plausible theories for IMBH formation in glob- 
ular clusters, or to assess whether globular clusters may 
follow the same black hole demographics correlations as 
galaxies. 

Isotropic models are not accurate for to Centauri, 
since they do not reproduce the observed deviations of 
[£ pm t/£ pmr ](i?) from unity. Use of isotropic models to 
fit a(R) yields spuriously high IMBH masses. Isotropic 
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models do provide a reasonable guide for calculation 
of predicted velocity distribution shapes (as opposed to 
RMS values) in the central region. The predicted Gauss- 
Hermite moments as function of radius for such mod- 
els fit the observations in the central arcmin to within 
\Ahi\ < 0.01. In models with an IMBH, the wings of the 
velocity distribution for an aperture on the very center 
are predicted to be more extended than in models with- 
out an IMBH. However, for cusp models with plausible 
IMBH masses we would not have expected more than a 
single star in the data set to have a velocity component 
in excess of 60 kms" 1 (~ 3.5 times the cluster velocity 
dispersion). No such stars were observed in our data. 
Given the small number statistics, this does not provide 
meaningful new constraints on the IMBH mass for this 
particular cluster. 

The dark mass invoked in cusp models may in principle 
be extended. The \ 2 °f the fit yields a limit adark < 
7" = 0.16 pc on the extent of any dark cluster at la 
confidence. The corresponding central density is > 500 
times larger than the central density of luminous matter. 
This would be indicative of quite an extreme amount of 
mass segregation. If u> Centauri does have a cusp and 
a dark mass, then an IMBH may be easier to explain 
theoretically than such a concentrated dark cluster. 

The IMBH mass suggested by NGB08 is strongly ruled 
out in all models. They found M BH = 4.0l% 5 x 10 4 M Q 
from isotropic modeling and Mbh = (3.0 ± 1.0) x 10 4 M Q 
from anisotropic modeling. By contrast, we find that 
any IMBH in excess of 1.8 x 10 4 M© is ruled out at 3a 
or higher confidence. This difference is not surprising, 
given that our observations have not reproduced the ar- 
guments on which their IMBH mass was based (namely, 
their choice of cluster center, their finding of an increase 
in dispersion towards that center, and their finding of a 
cusp around that center). If one ignores these fundamen- 
tal differences, then maybe our results do not appear so 
different from those NGB08. When comparing similar 
models between the two studies (i.e., cusp models with 
and without anisotropy), the differences in the implied 
IMBH masses are at the ~ 2a confidence level, given the 
uncertainties quoted by NGB08 (with our masses being 
lower). However, the more important difference is that 
we find that core models with no IMBH are perfectly 



consistent with the data. By contrast, NGB08 concluded 
that models with a core and models without an IMBH 
were ruled out with least 3a confidence. So while NGB08 
argued for a significant IMBH detection, we find by con- 
trast that there is no need to invoke an IMBH in Omega 
Cen at all. Our finding that uj Centauri cannot have 
an IMBH as massive as suggested by NGB08 is consis- 
tent with arguments presented by Maccaronc & Servil- 
lat (2008). They found that accretion models with the 
NGB08 IMBH mass predict much more X-ray and radio 
emission than is observed from the cluster center. 

To further strengthen the conclusions from this paper 
and our understanding of lo Centauri, it will be use- 
ful to proceed with the construction of more sophisti- 
cated equilibrium models. Methods for numerical cre- 
ation of axisymmctric models through orbit superposi- 
tion have been developed by, e.g., vdV06, van den Bosch 
et al. (2006), and Chaname et al. (2008). These mod- 
els can fit the observed differences between major and 
minor axis kinematics (discussed in Appendix B). They 
can also make use of the 47,781 high-quality proper mo- 
tions along the major axis at R w 1-6 arcmin, which 
were excluded from consideration here (discussed in Ap- 
pendix C). This will yield constraints on the inclination, 
rotation properties, and the presence of any kinematic 
subcomponents. It should also yield further improved 
limits on the possible presence of any dark mass in the 
center. Construction of non-equilibrium models for u> 
Centauri would be useful as well. A-body, Monte-Carlo, 
or Fokker-Planck models might shed light on the evo- 
lutionary status of lo Centauri, and on how it arrived 
at its present structure. This would also provide a more 
self-consistent description of the joint effects of mass seg- 
regation and energy equipartition than we have provided 
here. Moreover, it would be able to explicitly account 
for the Brownian motion of a possible IMBH. 



We are grateful to Ivan King, Dean McLaughlin, Pat 
Seitzer, and Glenn van de Ven for useful discussions 
and/or providing data in electronic format. Suggestions 
from the referees helped us improve the presentation of 
our results. 



APPENDIX 

A. EXTRACTING KINEMATICS FROM DISCRETE DATA 

To extract kinematics from a set of discrete data points, we select the stars that fall in a given area on the sky 
(e.g., a circular aperture, annulus, polar wedge, or square aperture). For these stars we then select the measurements 
Wi ± Awi in a given coordinate direction (e.g., line of sight, proper motion radial, proper motion tangential, proper 
motion major, or proper motion minor). The Wi may be either velocities in km/s or proper motions in mas/yr. The 
goal is to determine the underlying distribution from which the measurements are drawn. We assume that this is 
a Gaussian with mean W and dispersion a, and we set out to determine W and a with their error bars. If the 
measurement errors Awi are all identical then the problem has straightforward analytical solutions. However, this is 
not generally the case, so that more complicated analysis is required (e.g., Appendix A of vdV06). 

We start with a trial estimate of a. Each measurement Wi is then drawn from a Gaussian with dispersion di = 
(Awf + cr 2 ) 1 / 2 . Hence, the velocity W and its error bar AW follow from straightforward weighted averaging 

S! = ]TiM 2 , s 2 = ^2 w i/ d l w = s 2 /s 1 , aw = s; 1/2 . (Ai) 

i i 

We then determine the value of a which maximizes the likelihood 

L(a) = H^na 2 )- 1 / 2 exp{-K - W) 2 /(2[Aw 2 + a 2 ])}, (A2) 



2G 



van der Marel & Anderson 



by numerically solving the equation dL/da = 0. Having obtained ct, we return to an improved determination of W 
and we iterate the procedure until convergence. This provides a fixed point iteration scheme for determination of the 
joint maximum likelihood solutions for W and ct. 

As described in Appendix A of vdV06, the maximum likelihood solution for a yields a subtly biased estimate of 
the true dispersion. vdV06 used an analytical approximation to correct for this. Here we have used a Monte-Carlo 
procedure instead. After the determination of (W,a), we draw pseudo-data sets from the corresponding Gaussian 
probability distribution. To this we add random Gaussian errors drawn from the observational uncertainties. We 
then analyze numerous such pseudo-datasets in the same fashion as the real data. The statistics of the Monte-Carlo 
results provide both an estimate of the bias in a (which we use to correct our maximum likelihood estimate) and of 
the uncertainty Act. 

After the calculation of W and a has been completed, the RMS value a and its uncertainty can be obtained from 

a = {W 2 + a 2 ) 1/2 , Act = [{W AW) 2 + (ct Aa) 2 } 1/2 /a. (A3) 

For the HST proper motion data, our analysis procedure removed all mean motions (see Section 15.2. 3|) . Therefore, we 
know that the mean of the underlying distribution for any data set extracted from it should be W = 0. For this case 
we determined a more efficiently as the solution of the maximum likelihood equation dL/da = at fixed W = 0. 

B. MAJOR VERSUS MINOR AXIS 

In a spherical system one expects the projected kinematics to be the same along any axis through the center. One 
also expects that the projected kinematics for the proper motion coordinate along any fixed direction, when averaged 
over a circle, will be independent of the chosen direction. These expectations are not realized in practice, and this 
provides a kinematical signature of the axisymmetry of u> Centauri. vdV06 discussed this in detail for the ground-based 
kincmatical data. Here we briefly discuss this issue for the HST proper motion data. 

We have extracted major and minor axis subsamples from our HST proper motion sample, consisting of the stars 
lying within 10° from either axis. When averaged over all radii R <_71.7", we find for the radial proper motion 
direction that £ pmr = 0.758±0.010masyr _1 along the major axis, and £ pmr = 0.739±0.010masyr~ 1 along the minor 
axis. For the tangential proper motion direction we find that £ pm t = 0.691 ± 0.009 mas yr -1 along the major axis, 
and Sp m t = 0.764 ± 0.010 mas yr -1 along the minor axis. In a spherical system the kinematics would be the same 
along the two axes. In reality, E pmr is (2.5 ± 1.8)% lower along the minor axis than along the major axis, while S pm t 
is (10.6 ± 2.1)% higher along the minor axis than along the major axis. There is no statistically significant variation 
with radius in these comparisons (for R < 71.7"). 

An alternative way to look at these results is to consider the fixed major and minor directions on the sky, and to 
measure the RMS projected proper motion components S pmm j and S pmm ; along these directions, respectively. For 
the major axis subsample we find S pmm j/S pmm i ~ S pmr /Spurt = 1.097 ± 0.021. For the minor axis subsample we find 
Spmmj/Spmmi « S pmt /S pmr = 1.039 ± 0.020. When averaged over the whole sample (i.e., over all position angles) we 
find that S pmm j/fi pmm i = 1.048 ± 0.007. Therefore, the RMS projected velocity along the major axis is always larger 
than that along the minor axis. This was also found by vdV06, who obtained for the 718 stars with all three velocity 
components measured that £ pmm j/E pmm i = 1.055 ± 0.053 (see their Appendix C and also their Figure C.l). 

These results are relatively easily understood heuristically. For an axisymmetric_system that is not too far from 
edge-on, £ pmm j is primarily related to the pressure along the equatorial plane and £ pm mi is primarily related to the 
pressure along the symmetry axis. The former must exceed the latter in an oblate system, as dictated by the tensor 
virial theorem (e.g., Binney & Tremaine 1987). 

C. PROPER MOTION KINEMATICS FOR THE MAJOR AXIS HST FIELD 

As discussed in Section 15.2.11 there is also proper motion data available from Paper I for another field that was 
observed at multiple epochs with HST. This field lies roughly along the major axis, and therefore does not provide 
coverage of the full range of position angles in u) Centauri. We therefore did not use it for our model fitting. However, 
we address here the major axis profiles that can be extracted from these data. This allows a comparison with the 
central field and our model predictions to check for consistency. 

We defined both for the central field and the major axis field a major-axis subsample consisting of stars lying within 
10° from the major axis. This is the largest wedge for which coverage over the full range of allowed position angles is 
guaranteed at the majority of radii R. For the central field we have complete coverage for R < 115". For the major 
axis field we have complete coverage for 108" < R < 292". We extracted the profiles of E pmr (i?) and E pmt (i?) for 
both subsamples, following the procedures described in Section f5. 2. 21 The results are shown in Figure [TT1 The HST 
profiles are continuous at the boundary between the central and major axis fields, as they should be. 

At radii R < 71.7", we have complete coverage in the central field over all position angles. We can therefore calculate 
the ratios of the projected RMS proper motions measured within 10° from the major axis and averaged over an entire 
circle, respectively. This yields (£ pmr ) maj /(I] pmr ) C i rc = 1.011 ± 0.014 and (S pm t)maj/(S pmt )circ = 0.938 ± 0.014. If 
we assume as a simple approximation that these ratios are independent of radius, then this allows us to compare 
the predictions of our dynamical models to the observed major axis profiles. The curves in Figure [TT1 are the model 
predictions for our best-fitting anisotropic model with a core and no dark mass (see FigureJT^i), scaled by the previously 
quoted ratios (S pmr ) ma j/(S pmr ) c i rc and (£ pm t)maj/(E pm t) c irc- 
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Fig. 11. — Observed major axis profiles of the RMS projected proper motions E pmr and £ pm t in mas/yr, as function of projected radius 
R from the cluster center. Symbols are the same as in Figure |3p. For the HST data points we included stars within 10° from the major 
axis. For the vdV06 data points (magenta in the on-line version, with horizontal error bars) we included stars within a slightly larger angle 
from the major axis (17.5°), to increase the sample size. The vdV06 motions were multiplied by a factor 1.023 to correct to the same 
characteristic main-sequence stellar mass as for the HST sample (see Section 15.1.31 . HST data points at R < 110" (closed dots) are based 
on the data from the central field, and those at R > 110" (open dots) are based on the major axis field. The figure focuses on the central 
region and uses a linear scale in R, by contrast to Figure [3p. The curves are predictions based on the best-fitting anisotropic model with 
a core and no dark mass, as discussed in the text. The data for the major axis field are consistent with expectation, given the detailed 
modeling of the central field and the inability of the HST data to measure absolute rotation in the pmt direction. 

The model curves thus obtained provide an excellent fit to the S pmr (i?) data for the major axis field. Also, the 
Sp m r(-R) data for the major axis field are consistent with the data from vdV06, which are also shown in the figure. By 
contrast, the fit to the £ pmt (i?) data for the major axis field is good only for R < 200". At larger radii, the data fall 
below the model, and they also fall below the data from vdV06. We attribute this to the inability of the HST proper 
motion data to measure absolute rotation V pm t. As a result, we expect to progressively underestimate £ pm t with 
increasing radius, since V pm t/S pm t increases with radius (see Section [5.2. 3p . The fact that the data-model comparison 
is more successful for E pmr (i?) is explained by the fact that pmr is orthogonal to direction of absolute rotation, and it 
is therefore insensitive to it. 

Taking into consideration the limitations inherent to the data for the major axis field, and the results presented 
here, we have no reason to believe that detailed modeling of the data for the major axis field would significantly alter 
the conclusions of our paper. 
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